DROMPA3: その1 インストール

今回からは私が開発したDROMPA3の利用法について解説します。

DROMPAとは

ChIP-seq解析のためのパイプラインツールです。ピーク抽出の他、品質評価、可視化、複数サンプルの比較解析などができます。 複数のサンプルを同時に解析できること、pdf形式でデータを出力できること、ChIP-seq解析の様々なステップをオールインワンで含んでいることが特長です。 以下が元論文です。

DROMPA: easy-to-handle peak calling and visualization software for the computational analysis and validation of ChIP-seq data - Nakato - 2013 - Genes to Cells - Wiley Online Library

この論文はversion 1の頃のもので、現在のversion3は色々な部分が変わっているので、詳細はDROMPA3付属のマニュアルを読んでいただくと良いと思います。
ちなみに名前の由来は、DRaw and Observe Multiple enrichment Profiles and Annotation の頭文字を取ったもの…ということになっていますが、 本当はFC東京のマスコットのドロンパから名前をもらいました。（FC東京サポなので）

以下、インストール方法です。

関連ライブラリのインストール

まずインストールに必要となるライブラリ群をインストールします。詳細はDROMPA3のWebサイトを参照してください。

Ubuntuの場合:

 $ sudo apt install git gcc libgtk2.0-dev libgsl-dev

CentOSの場合:

 $ sudo yum -y install zlib-devel gsl-devel gtk2-devel

ダウンロード・コンパイル

githubからダウンロードし、以下のようにしてコンパイルします。

 $ git clone https://github.com/rnakato/DROMPA3
 $ cd DROMPA3
 $ make

関連ツールのインストール

BAMフォーマットのファイルを入力にする場合、SAMtools が必要になります。 インストール方法は以下の記事を参照してください。

SAMtoolsとリダイレクト - Palmsonntagmorgen

また、複数のpdfを結合するために Coherent PDF というツールを使いますので、こちらもgithubからダウンロードします。

 $ git clone https://github.com/coherentgraphics/cpdf-binaries

PATHを通す

最後にPATHを通します。PATHとは何ぞや？という方は以前の記事を参照してください。
例としてここではホームディレクトリに my_chipseq_exp というディレクトリを作り、その中にDROMPA3とCoherent PDFをダウンロードした場合のPATHを記載します。

 $ export PATH = $PATH:$HOME/my_chipseq_exp/DROMPA3:$HOME/my_chipseq_exp/cpdf-binaries/Linux-Intel-64bit/

これでインストールは完了です。 以下のコマンドをタイプして、DROMPA3のヘルプが表示されれば成功です。

 $ drompa_draw
Usage: drompa_draw [--version] <command>
Command: PC_SHARP    peak-calling (for sharp mode)
         PC_BROAD    peak-calling (for broad mode)
         PC_ENRICH   peak-calling (enrichment ratio)
         GV          global-view visualization
         PD          peak density
         FRIP        accumulate read counts in bed regions specified
         CI          compare peak-intensity between two samples
         CG          output ChIP-reads in each gene body
         GOVERLOOK   genome-wide overlook of peak positions
         PROFILE     make R script of averaged read density
         HEATMAP     make heatmap of multiple samples
         TR          calculate the travelling ratio (pausing index) for each gene

SAMtoolsワンライナー覚書

順次追加するかも。versionは1.5です。
.sort.bam はソート済BAMを表します。

SAM -> BAM 変換

 $ samtools view -bS sample.sam > sample.bam

BAM -> SAM 変換

 $ samtools view -h sample.bam > sample.sam

BAMをソート

 $ samtools sort sample.bam > sample.sort.bam

SAM -> ソート済BAM

 $ samtools view -bS sample.sam | samtools sort > sample.sort.bam
 または
 $ samtools sort sample.sam > sample.sort.bam

BAM indexを作成

 $ samtools index sample.sort.bam

複数のBAMファイルをマージ

 $ samtools merge sample.merged.bam sample1.bam sample2.bam

SAM/BAMに含まれるリード数をカウント（unmappedも含む）

 $ samtools view -c sample.bam

SAM/BAMに含まれるリード数をカウント（mapped readのみ）

 $ samtools view -F 0x04 -c sample.bam

マップされなかったリード数をカウント

 $ samtools view -f4 -c sample.bam

染色体別のマップリード数をカウント（index要）

 $ samtools idxstats sample.bam

もうちょっと詳細な情報が知りたい

 $ samtools flagstat sample.bam

さらに詳細な情報が知りたい

 $ samtools stats sample.bam

BAMからmapped readのみを抽出

 $ samtools view -F 0x04 -b sample.bam > sample.mapped.bam

BAMからPCR duplicateを除去

 $ # single-end
 $ samtools rmdup -s sample.sort.bam sample.sort.rmdup.bam  
 $ # paired-end
 $ samtools rmdup sample.sort.bam sample.sort.rmdup.bam    

BAMからPCR duplicate除去済のmapped readを抽出

 $ samtools view -F 0x04 -b sample.sort.bam | samtools rmdup -s - sample.rmdup.sort.bam

BAMから１番染色体にマップされたリードのみを抽出しfastqとして出力

 $ samtools view -b sample.bam chr1 | samtools fastq

viewer (tview)

 $ samtools tview aln.bam ref.fasta  # using reference genome
 $ samtools tview aln.bam            # not using reference genome

• .: forward, identical to reference genome
• ,: reverse, identical to reference genome
• ACGTN: forward, mismatch to reference genome
• acgtn: reverse, mismatch to reference genome

pileup format

 $ samtools pileup aln.bam -f ref.fa   # using reference genome
 $ samtools pileup aln.bam             # not using reference genome

tviewと違い、何もマップされない領域は表示されない。

環境変数PATHの通し方

同内容の記事はたくさんありますが、やはり避けては通れないので…

環境変数PATHとは

githubなどからツールを新たにダウンロードした場合、その実行ファイルを起動するには実行ファイルのありかを直接指定する必要があります。
$ ./bowtie2-2.2.9/bowtie2 のような感じですね。
しかし ls や pwd のようにLinuxにもともと含まれるコマンドや、ssh や R のようにapt-getなどの方法でインストールしたコマンドは、パスを指定する必要なく、 ls ssh とタイプするだけで起動します。

なぜでしょうか。

結論から言うと、Linuxには環境変数PATHというものがあり、PATHに指定されているディレクトリ内の実行ファイルは自動的に検索されるため、絶対パスで指定しなくてもよいのです。

実行ファイルの所在を調べる

では、ls や R の実行ファイルはどこにあるのでしょうか。 whichコマンドか whereコマンドを使うと調べることができます。

$ which ls
/bin/ls
$ which R
/usr/bin/R

ls は/bin, R は/usr/bin 内にあることがわかりました。 Linuxに元々付属するコマンドは/bin, aptなどでインストールしたものは/usr/binに含まれる傾向があります。 なお、特定のユーザのみのためにインストールされるツールの場合は /usr/local/bin にインストールされることもあります。

環境変数PATHを表示する

echo $PATH とタイプすると、環境変数PATHに登録されているディレクトリを一覧表示できます。

$ echo $PATH
/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:

/bin や /usr/bin が含まれていることが確認できます。 ここに表示されていないパスについては、実行ファイルを直接指定する必要があるということですね。

PATHを一時的に追加する

環境変数PATHに新たに参照ディレクトリを追加したい場合（PATHを通す と言います）、 一時的に追加する場合はコマンドライン上で以下のようにタイプします。
（ここでは上述のbowtie2を例とします）

$ export PATH=$PATH:/home/rnakato/bowtie2-2.2.9/

こうすると/home/rnakato/bowtie2-2.2.9/にPATHが通り、bowtie2とタイプするだけで実行できるようになります。 複数PATHを通したい時にはコロン(:) で挟んで併記します。

$ export PATH=$PATH:<パス1>:<パス2>:<パス3>:...

ちなみに$PATHと最初に記載しているのは、既存のPATHを上書きしてしまわないためです。
この方法でPATHを通した場合、そのターミナルを閉じると追加したPATHも消えてしまいます。

PATHを永続的に追加する(.bashrc)

常にPATHを通したい場合はホームディレクトリにある .bashrc (または.bash_profile)に以下の行を追記します。 （ファイルが存在しない場合は新規作成で問題ありません）

export PATH=$PATH:<パス1>:<パス2>:<パス3>:...

.bashrcに書かれた設定は次回ターミナル起動時から有効になりますので、既存のターミナルで新しく追加した設定を有効にしたい場合は以下のコマンドを実行します。

$ source ~/.bashrc

同一プログラムが複数箇所にある場合

ツールによっては、ソースからインストール、aptやyumでインストール、あるいはcondaのようなツールでインストールと、 複数の方法でインストール可能なものがありますが、インストールされるディレクトリはそれぞれ異なります。
複数のPATH上に同一プログラムがインストールされている場合、最初にヒットしたPATHが優先されます。
つまり、bowtie2が/usr/bin と /home/rnakato/bowtie2-2.2.9/ に存在し、PATHが /usr/bin:/home/rnakato/bowtie2-2.2.9/ となっていた場合、 単にbowtie2とタイプすると/usr/bin/bowtie2が起動します。 /home/rnakato/bowtie2-2.2.9/bowtie2 のように絶対パスで記述することで好きな方を起動可能ですが、バージョンが異なる同一ツールをあちこちにインストールするのはトラブルの元ですので、 プログラムは一箇所のみのインストールに留めるべきでしょう。
（一般にaptやyumでインストールされるプログラムはバージョンが古いことが多いです）
また、which コマンドを利用して、自分が起動しているプログラムがどこにインストールされたものかを常にチェックしておく癖をつけておくとより安心です。

Bowtie2が system CPU を大量に消費している件で

Bowtie2の最新バージョンは2017/10/10 時点でVersion 2.3.3.1となっています。

Change logを見ると、version 2.3.0 において major updateが施され、 マルチスレッドを使用した場合のスケーラビリティを改善したとあります。 具体的には、利用するライブラリをTBBに変更したようです。

http://bowtie-bio.sourceforge.net/bowtie2/index.shtml

ところがその後、bowtie2で -p 12 など多数のCPUを指定した場合に、Systemが多量のCPUを消費してしまい、 実際のマッピングの速度が非常に遅くなるという現象が報告されています。

bowtie2 use high sys %Cpu . · Issue #62 · BenLangmead/bowtie2 · GitHub

私の手元でも確認しています。

Systemがマッピングの倍くらいCPUを消費しています。

bowtie2のマニュアルによると、ソース（bowtie2-2.3.3.1-source.zip）をダウンロードしたうえで　NO_TBB=1オプションを付けて　makeでコンパイルすると TBBをオフにできるとのことなので、そのようにしたところ、この現象は改善しました。

major update前のVersion2.2.9を使っても改善しますが、古いバージョンはfastq.gz形式をそのままでは受け付けないので、上記の方法の方がよいと思います。 この現象は今年の2月時点で対応中ということになっていますが、現在の最新バージョンでも根本的解決には至っていないようです。

なお、version2.3.0以降で怪しげなバグが他にも報告されているので、まだ問題が残っているかもしれません。。。 しばらく様子を見たいと思います。

SAMtoolsとリダイレクト

SAMtools

先日紹介したリードのマッピングの記事では、出力をSAM形式に指定していました。

rnakato.hatenablog.jp

SAM形式はファイルサイズが非常に大きく読み込みにも時間がかかるので、バイナリ形式のSAMであるBAM形式が下流解析でよく利用されています。
ここではSAM形式をSAMtoolsを使ってBAM形式に変換します。

SAMtoolsのインストール

Ubuntuであればaptでインストールできます（ややバージョンが古い）。

 $ sudo apt install samtools

Anacondaを使っている人であればcondaでもインストールできます。

 $ conda install -c bioconda samtools

どちらでもない人はSAMtoolsのWebサイトからダウンロードしましょう。

BAM作成

以下のコマンドを実行すると、SAM形式のファイルを読み込んでBAM形式に出力します。

 $ samtools view -bS sample.sam > sample.bam

-b はBAM系形式で出力するオプションで、-Sは入力ファイルの形式を自動認識させるオプションです。-hをつけると、SAMファイルのヘッダを含めるようになります（デフォルトでは省略される）。
-@オプションをつけると使用CPU数を指定できます。

 $ samtools view -bS sample.sam -@ 4 > sample.bam

BAMのソート

上記コマンドで得られたbamファイルはソートされていませんが、多くのツールではsorted BAMを要求されます。
以下のようにしてbam をソートしましょう。

 $ samtools sort sample.bam > sample.sort.bam

これでソートされたbamファイルが生成されました。
ソート前のファイルは不必要なので削除しましょう。

注意点：古いバージョンのsamtools (0.1.19など)とはsortコマンドの指定方法が変わっていますので、 古いバージョンを想定したスクリプトを新しいバージョンのSAMtoolsで動かそうとするとエラーになります。 sortコマンドがエラーになる時はバージョンを意識してみましょう。

リダイレクトによるsorted bam生成

上記のコマンドはBAM作成 -> sorted BAM作成を別々に行っており、中間ファイルを削除する手間もあるのが面倒です。 そこで、以下のようにリダイレクトで生成されるBAMファイルを直接ソートしてみましょう。

 $ samtools view -bS sample.sam -@ 4 | samtools sort > sample.sort.bam

| は、左側のコマンドから標準出力に出力される結果をそのまま右側のコマンドの入力にするという命令です。 samtools viewの出力を直接samtools sortに投げることで、中間ファイルを作成することなくsorted bamを作成することができます。

なお、以下のようにsortコマンドに直接SAMファイルを指定しても同じことができるようです。

 $ samtools sort -@ 4 sample.sam > sample.sort.bam

注意点：samtoolsは内部で生成される中間ファイルをカレントディレクトリに生成しますので、 一度の多量のファイルを並列で変換すると、中間ファイルが衝突して結果がおかしくなる可能性があります。

BAMのindex作成

一部のツールはBAMの読み込みを高速化するためにBAMのindex（索引）ファイルを要求します。
以下のコマンドでindexを作成します。なお、indexコマンドはsort済のBAMファイルしか受け付けません。

 $ samtools index sample.sort.bam -@ 4

こうすると、sample.sort.bam.bai　という名前のindexファイルが生成されます。

BAMファイルを作成する時は、sortとindex作成を常に行うようにしておけば間違いないでしょう。

おまけ：マッピングツールの出力をそのままsorted bamに変換する

マッピングツールの出力をリダイレクトでsamtoolsに投げてやると、そもそもSAMファイルを作成・削除する手間が要りません。

 $ bowtie2 -x index sample.fastq -p4 | samtools view -bS - | samtools sort > sample.sort.bam

for文でたくさんのサンプルを処理する時も１行で書けて楽ちんです。

 $ for num in $(seq 1 20); do
 $    bowtie2 -x index sample$num.fastq -p4 | samtools view -bS - | samtools sort > sample$num.sort.bam
 $ done

バイオインフォマティクスのためのpython環境構築方法を考える (10/13追記あり)

先日、以下の記事でLinux上でのpython環境構築にはpyenvが良いと書いたのですが、

pyenvでPython環境を構築する - Palmsonntagmorgen

少し気が変わってきました。

以下、現状の考えをまとめます。
私もPythonにはそこまで詳しくないので、良い方法をご存じの方がおられましたらご教示ください。

目的

Python系の既存ツールを使えるための環境づくり。
既存ツールには2系と3系で書かれたものがあるので、どれでもインストールできるようにするためには、 2系と3系の環境を共存させる必要がある。

（ので、例えば2系のツールしか当面必要ありませんという人ならそもそもpyenvが要りません）

現在のツールの状況

今のところ、NGS周辺の解析ツールは2系で書かれたツールが多数を占めています。
たとえばMACS2は2系でしか動きません。
2系のサポートは2020年までなのでいずれは3系に移行していくのでしょうが、 現時点ではどちらか一つ選ぶとしたら2系にせざるを得ないでしょう。
（開発者のやる気を考えると、現在2系で書かれているツール群の大半は3系に書き直されることはないでしょうし…）

一方、Anacondaを想定したツールも増えてきており(biocondaとか)、これからPythonで開発される解析ツールはAnacondaベースになっていくのではないかと思います。
なので、やるべきはAnaconda2とAnaconda3の共存ということになります。

問題1

pyenvを使う使わないに関わらず、Anacondaをインストールした際にはローカルのpythonより優先させるために 環境変数PATHの上位にanacondaを記載すると思います。 そうすると、以下の記事にもあるように、Anacondaが持っているツールがOSの持っているツールを覆い隠してしまいます。

pyenvが必要かどうかフローチャート

つまり例えばAnacondaでRがインストールされると、以降はRを起動した時にそちらが呼び出され、ローカルのRは使われなくなってしまいます。

問題2

Anacondaは今のところそれほど安定ではなく、ローカルへのインストールでは起きなかったエラーが起きることがよくあります。
2017/9/13現在、最新のAnaconda2-4.4.0を使ってRをインストールすると、（環境によると思いますが） install.packagesで以下のようなエラーが出るようです。（Rcppの例）

2017/10/13 追記：　下記エラーはgccのバージョンを4.8にダウングレードすることで解消されました。

** testing if installed package can be loaded
Error: package or namespace load failed for ‘Rcpp’ in dyn.load(file, DLLpath = DLLpath, ...):
 unable to load shared object '~/.pyenv/versions/anaconda2-4.4.0/lib/R/library/Rcpp/libs/Rcpp.so':
  ~/.pyenv/versions/anaconda2-4.4.0/lib/R/library/Rcpp/libs/Rcpp.so: undefined symbol: _ZNSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEED1Ev
Error: loading failed
Execution halted
ERROR: loading failed

Rcppだけならばcondaで対処可能ですが、setup.py install でツールをインストールした場合でも類似の現象が起きる場合があります。
ローカルとAnaconda内のコンパイラ・ライブラリのバージョンの相違に起因するようですが、 ローカルのRならば考えなくて良いエラーの対処に多くの時間を取られてしまうのはかなりのストレスです。

問題3

ならばAnacondaでRをインストールしなければよい、という話なのですが、解析ツールのいくらかはdependenciesにRを持っており、 そのツールをインストールすると自動的にRもインストールされてしまいます。(MAGeCK-VISPRなど)
fastx-toolkitをインストールすると何とperlがインストールされます。
（perlに関してはplenvを使うことで対処可能です）

それらのツールのインストールにはAnacondaを使わないで…と回避することは現状可能ですが、 今後もAnacondaに依存したツールが増えるであろうことを考えると、どうするか検討が必要です。

問題4

共有サーバ上でpyenvを用いて統一環境を構築した場合の問題ですが、
他の一般ユーザはpyenvのlocalや仮想環境で2系と3系を切り替える権限を持たないため、 ある人は3系のツールを使いたい、ある人は2系のツールを使いたい…という風に使い分けるのがかなり困難です。
（本来はユーザごとに自力でpyenvをインストールして環境設定すべきなのですが、それを要求するのは難しいでしょう）
また、pyenvを使ってAnaconda2と3を共存させる限り、~/.pyenv/versions/anaconda*/bin/ にパスを通すことができません。
(どちらか片方だけを使うつもりであるなら通してもよいのですが、そうすると共存とは言えないし…)

さらに、pyenvを利用して2系と3系を共存させ、 かつ3系ではRがインストールされているが2系ではインストールされていないというような場合、 2系を利用している時にRを起動しようとすると「Rがインストールされていません　3系の環境なら利用可能です」とエラーが表示され、 ローカルのRを呼び出してくれません。

そもそも、1カ所で統一できるはずのRを何箇所にもインストールすること自体に違和感があります。

対処

以上のことを考えると、サーバ上でpython共有環境を構築する際には、 今のところ以下の方法が良いのではないかと思っています。

• 生のAnaconda2をインストールし、パスを通す。
• AnacondaではRをインストールせず、ローカルのR (/usr/bin/R)に統一。
• 別にAnaconda3もインストールするが、パスは通さず、必要な時だけ絶対パスで起動。

このようにした場合、pythonでは基本的にAnaconda2が呼び出され、RはローカルのRが呼び出されます。
なおpyenvを使ってAnacondaをインストールし、localや仮想環境を使わない場合でも同様の環境は構築できますが、 生Anacondaの方がどちらを利用しているかはっきりするのが良いかと思いました。

Rに依存しているAnacondaツ―ルがローカルのRを正しく呼び出してくれるのかは未確認です。。

もっとよい方法があれば教えてください。。

雑感

これからはPython3系になっていくので3系のみインストールすれば十分、という記事はよくありますが、 2系で書かれたツールを使わざるを得ないバイオインフォの世界ではそうもいかないのが辛いところです。

自分が利用したいツールが何かによって最適な環境はかなり変わるはずなので、 自分の条件に沿うやり方を探すのがよいと思います。
私はかなり色んなツールを要求されるのでこんなことになっていますが…

Anaconda内のRやらのインストールでエラーが出なければ、特に問題はないんですけどね。。早く安定になって欲しい。

Readをゲノムにマッピング (その3) 圧縮ファイルを入力にする方法

マッピングの記事その３。

Readをゲノムにマッピング (その１) - Palmsonntagmorgen
Readをゲノムにマッピング (その２) - Palmsonntagmorgen

圧縮ファイル(fastq.gz)を直接マッピングの入力にする方法です。

圧縮ファイルのままマッピングしたい

fastqファイルは容量が大きいので、数が増えてくるとHDDを圧迫します。
なので.gzで圧縮した状態で保管しておきたいですね。 しかし圧縮されたfastq.gzファイルをbowtieにかけるとエラーで終了してします。

 $ bowtie -S hg38 sample.fastq.gz -n2 -m1 -p4 > sample.hg38.sam
Error: reads file does not look like a FASTQ file
terminate called after throwing an instance of 'int'
Aborted

これはbowtieが圧縮ファイルの入力に対応していないからですが、 以下のようにコマンドを修正すると、マッピングが通るようになります。

 $ bowtie -S hg38 <(zcat sample.fastq.gz) -n2 -m1 -p4 > sample.hg38.sam

これはコマンドの中でgzファイルをzcatで展開し、その出力をbowtieの入力として投げる、というコマンドです。
bowtieに限らず、圧縮ファイルを受け付けないツール類（の一部）は、この方法で.gzファイルを直接入力にすることが可能です。

ちなみにbowtie2 (version 2.3.1 以降), bwaはそのままでもマッピングが通ります。

 $ bowtie2 -x hg38 sample.fastq.gz -p4 > sample.hg38.sam

 $ bwa mem -t 4 hg38 sample.fastq.gz > sample.hg38.sam