Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

記事一覧

解析環境構築

環境変数PATHの通し方

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

バイオインフォマティクスのためのpython環境構築方法を考える

データ生成

常染色体と性染色体のみのゲノム配列ファイル genome.fa を作成する - Palmsonntagmorgen

2bit genome を作成する - Palmsonntagmorgen

genome tableを作成する - Palmsonntagmorgen

 

fastqデータ取得

SRAからfastqを取得する - Palmsonntagmorgen

ENA,DDBJからfastqを取得する - Palmsonntagmorgen

 

ゲノムマッピング

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

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

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

 

環境変数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

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

f:id:rnakato:20171010172635p:plain

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

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

f:id:rnakato:20171010173010p:plain

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文でたくさんのサンプルを処理する時も1行で書けて楽ちんです。

 $ 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が必要かどうかフローチャート - Qiita

つまり例えば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) 圧縮ファイルを入力にする方法

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

Readをゲノムにマッピング (その1) - Palmsonntagmorgen
Readをゲノムにマッピング (その2) - 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

Error: libpng12.so.0 not found (VMware)

メモ。

VMwareにインストールしたUbuntu上で ftp://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/UCSC toolsを実行すると以下のエラーが出る。

 error while loading shared libraries: libpng12.so.0: cannot open shared object file: No such file or directory

libpng12.so.0 は VMwareでは(手元の環境では)/usr/lib/vmware-tools/lib/libpng12.so.0/ にあるので、以下のコマンドでリンクを張ると動作するようになる。

 $ sudo ln -s /usr/lib/vmware-tools/lib/libpng12.so.0/libpng12.so.0 /usr/lib/x86_64-linux-gnu/libpng12.so.0

参考:

apt - Where can I find package 'libpng12.so.0'? - Ask Ubuntu