Palmsonntagmorgen

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

記事一覧

解析環境構築

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

 

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

先日、以下の記事で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の例)

** 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, 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

R in Anaconda にrJavaをインストール(8/15追記)

個人的につまったので忘備録。

anaconda は Linux (Ubuntu 16.04)上でpyenvを使ってインストールしているという前提です。

 $ conda create -c r r-irkernel # anaconda内にRのインストール

Rを起動してrJavaをインストールしようとするとエラーになる。

 > install.packages("rJava")

 ...
 gcc -std=gnu99 -o libjri.so Rengine.o jri.o Rcallbacks.o Rinit.o globals.o rjava.o  -shared -L/usr/lib/jvm/java-8-openjdk-amd64/jre/lib/amd64/server -ljvm -Wl,--export-dynamic -fopenmp  -L/home/rnakato/.pyenv/versions/anaconda3-4.4.0/lib/R/lib -lR -lpcre -llzma -lbz2 -lz -lrt -ldl -lm -liconv -licuuc -licui18n
 /usr/bin/ld: -liconv が見つかりません
 collect2: error: ld returned 1 exit status
 Makefile.all:38: ターゲット 'libjri.so' のレシピで失敗しました

-liconv を認識できずに色々苦労していたのだが、実はrJavaもcondaでインストールできることに気づいた。

 $ conda install -c r r-rjava

だがRを起動してrJavaをロードすると libjvm.so がないという理由でエラーになる。

 > library("rJava")
 Error :  .onLoad は loadNamespace()('rJava' に対する)の中で失敗しました、詳細は:
  call: dyn.load(file, DLLpath = DLLpath, ...)
  error:  共有ライブラリ '/home/rnakato/.pyenv/versions/anaconda3-4.4.0/lib/R/library/rJava/libs/rJava.so' を読み込めません:
  libjvm.so: 共有オブジェクトファイルを開けません: そのようなファイルやディレクトリはありません
 エラー:  ‘rJava’ に対するパッケージもしくは名前空間のロードが失敗しました

手元の環境ではlibjvm.soは/usr/lib/jvm/java-1.8.0-openjdk-amd64/jre/lib/amd64/server にあるので、LD_LIBRARY_PATH を通す。

 export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/lib/jvm/java-1.8.0-openjdk-amd64/jre/lib/amd64/server 

これでrJavaをロードできるようになりました。

なお、R CMD javareconf を使う方法は手元の環境ではうまくいきませんでした。

8/15 追記

conda install -c r r-rjava でrJavaをインストールしても、 R内でinstall.package(“rJava”)すると再インストールが始まり、同様のエラーとなる。
そのためrJavaに依存しているxlsなどをインストールできないが、 library(“rJava”)でロードしてからインストールするとエラーなくインストールが完了する。

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

前回の続きです。

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

今回ではbowtie, bowtie2, bwaのマッピングコマンドを説明します。

どのマッピングツールも、ゲノム配列をindex配列にまず変換し、そのindexに対してマッピングするという手順を踏みます。 直接ゲノム配列を検索するのではなく、indexを用いることで、高速・省メモリに計算することができます。
現在は、Burrows-Wheeler変換 (Burrows-Wheeler Transform; BWT) を用いてゲノムをindex化する手法が主流になっており、 上記3ツールは全てこの手法に基づいています。

Bowtieでマッピング

Bowtie, Bowtie2のindexはWebからダウンロードすることもできますが、ここではゲノム配列から自分で作成することにします。 human genome build hg38 を例とします。
index作成時のツールのバージョンとマッピング時のバージョンが異なる場合、indexの互換性が失われる可能性がありますので、 どのバージョンでindexを作成したかは記録しておきましょう。

index作成

 $ bowtie-build --version  # バージョン確認
 $ bowtie-build hg38.fa hg38 # index作成
 $ bowtie-build -C hg38.fa hg38-cs # index作成 (colorspace用)

colorspace dataを扱うことがなければ、3行目は必要ありません。

マッピング

 $ bowtie -S hg38 sample.fastq -n2 -m1 -p4 > sample.hg38.sam

-S で出力フォーマットをSAMに指定しています。 -p4 は CPUを4つ使うことを示します。 -n オプションは許容ミスマッチ数を示し、-n はBLAST(MAQ)風のアラインメントスコアで評価しますが、-v オプションを使うと単純にミスマッチ数のみで評価するようになります。
-m1 はユニークマッチのみ出力するオプションで、-m10 とすると、ゲノムの10箇所以下までマップされるリードを出力します。 -k1 オプションを入れると、1リードにつきベストスコアである1箇所のみを出力します。-n2 -m1の代わりに -n2 -k1とすると、マルチリードあり、1リードにつき出力は1箇所のみ、となります。 --best --strata オプションを使った場合もマルチプルリードも出力するようになり、出力されるのは1リードにつきベストスコアである1箇所のみになります。
colorspace dataをマッピングする場合は、-Cオプションをつけ、colorspace用のインデックスを指定します。

 $ bowtie -C -S hg38-cs sample.fastq -n2 -m1 -p4 > sample.hg38.sam

Bowtie2でマッピング

index作成

 $ bowtie2-build --version  # バージョン確認
 $ bowtie2-build hg38.fa hg38 # index作成

Bowtieとほぼ同じですね。

マッピング

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

indexは-x オプションで指定します。 Bowtie2はデフォルトの出力フォーマットがSAM形式のため、-Sオプションはありません。 出力ファイルにはマルチリードも含まれます。
マッピングのパラメータは色々あるのですが、それらの値をまとめた --very-fast, --fast, --sensitive, --very-sensitive というオプションが用意されています。デフォルトは--sensitiveです。
なおBowtie2は入力のreadにindelが入っている場合に分割して出力?するらしく、入力read数とSAMファイルのread数が異なる場合があります。

BWAでマッピング

index作成

 $ bwa  # バージョン確認
 $ bwa index -p hg38 hg38.fa # index作成

BWAは何故か--versionオプションがありません。。

マッピング

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

BWAもSAM形式で出力されます。出力ファイルにはマルチリードも含まれます。マッピングのパラメータは色々あるのですが、ちゃんと試したことはありません。。
PacBioとNanoporeのロングリード用のパラメータセッティングが用意されており、-xオプションで指定可能です。

その他

paired-end mappingについては記載していませんが、ほぼ同様の方法でマップすることができます。
基本的にどのツールを使っても問題ありませんが、RNA-seq, Hi-Cなどのツールの内部で特定のマッピングツールを要求されることもありますので、全て使えるようにしておくとより万全です。

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

NGS解析の最初のステップは、シーケンサから出力されたfastq形式のリード配列をゲノム配列にマップするマッピングです。 これにより、ゲノム上のどの領域から得られたリードなのかを知ることができます。

マッピングツール

ChIP-seq解析で用いられるマッピングツールは、bowtie, bowtie2, bwa の三種類です。 以下、それぞれの簡単な解説。

Bowtie

短いリード配列(〜50bp程度)を高速・高精度にマップします。 調べるのは配列ミスマッチのみで、挿入・欠失(indels)は考慮しませんが、 通常のマッピングには十分な精度を持ちます。 一方、長いリード配列(100bp以上など)にはあまり向いていません。
また、SOLiDシーケンサから出力されるcolorspace data のマッピングにこの中では唯一対応しています。 (バージョン1.2.0ではエラーが出たので、バグが修正されるまでは1.1.2を使うと良いでしょう)
colorspace dataは新たに生産されることはありませんが、過去のデータを再解析する際に必要になることがありますので、 その時はBowtieを使いましょう。
(というか、うちが以前メインで使っていただけなのですが。)

Bowtie2

bowtieの改良版で、挿入・欠失を考慮できるようになっている点が特徴です。 また、長いリード配列やペアエンド配列に対してbowtieよりも高精度です。
デフォルトでは リード配列全体をアラインメントするEnd-to-end alignment (グローバルアラインメント)を採用しますが、 --local オプションをつけると特に一致度の高い部分領域を探索するLocal alignmentに変わります。 localモードを試したことはないのですが、readとgenomeの配列一致度が低い場合には威力を発揮しそうです。
一方、colorspace read には対応していません。 また、unique read (ゲノムのただ一箇所のみにマップされるリード)と multiple read (ゲノムの複数箇所にマップされるリード)を分割することが難しいようです。(bowtieだと簡単)
Bowtie1と2のより詳細な違いが知りたい方は、Manualの"How is Bowtie 2 different from Bowtie 1?“を参照してください。

BWA

bowtie2と同じく挿入・欠失を考慮できます。 昔のバージョンではマッピングに aln -> samse/sampe と2段階踏む必要があり面倒だったのですが、 現在のバージョンではmemコマンドが実装され、より簡便になっています。
高精度ですが、こちらもunique read とmultiple read を区別することができず、常に両方出力されます。 また、マッピング後にスタッツを表示してくれない点が、Bowtieと比べるとやや不便です。
colorspace read に関しては、 csfastqを作成すればマッピング可能と書いてあったりするのですが、自分は成功したことがありません。 うまくいくやり方をご存知の方はご教示ください。。

どのツールを使えばいいか?

通常のChIP-seq解析においては、リード配列とゲノムの一致度が非常に高いので、 どのツールを使っても精度面ではほとんど差がありません。
ですので、やりたい作業(ユニークマッチだけ取りたい、挿入・欠失を考慮したいなど)によってツールを選ぶのがよいでしょう。
RNA-seqやPacBioなどのロングリードの場合はこの限りではありません。)
次回はマッピングのコマンドについて解説します。

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