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

前回の続きです。

Readをゲノムにマッピング (その１) (2017/12/19, 2018/11/19 追記あり) - Palmsonntagmorgen

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

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

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を扱うことがなければ、３行目は必要ありません。

マッピング

 $ 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 オプションを入れると、１リードにつきベストスコアである１箇所のみを出力します。-n2 -m1の代わりに　-n2 -k1とすると、マルチリードあり、１リードにつき出力は１箇所のみ、となります。 --best --strata オプションを使った場合もマルチプルリードも出力するようになり、出力されるのは１リードにつきベストスコアである１箇所のみになります。
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をゲノムにマッピング (その3) 圧縮ファイルを入力にする方法 - Palmsonntagmorgen

Readをゲノムにマッピング (その１) (2017/12/19, 2018/11/19 追記あり)

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

マッピングツール

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

Bowtie

短いリード配列（〜50bp程度）を高速・高精度にマップします。 調べるのは配列ミスマッチのみで、挿入・欠失（indels）は考慮しませんが、 通常のマッピングには十分な精度を持ちます。 一方、長いリード配列（100bp以上など）にはあまり向いていません。
また、SOLiDシーケンサから出力されるcolorspace data のマッピングにこの中では唯一対応しています。 （バージョン1.2.0ではエラーが出たので、バグが修正されるまでは1.1.2を使うと良いでしょう）

（2017/12/19追記：bowtieのgithubにリクエストを送ったところ、version1.2.2で上記エラーは解消されました。）

（2018/11/19追記：csfastqでなく、csfasta + qual に分かれたファイルをマップする時はやはりversion1.2.2だとエラーになるようです。version1.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 と２段階踏む必要があり面倒だったのですが、 現在のバージョンではmemコマンドが実装され、より簡便になっています。
高精度ですが、こちらもunique read とmultiple read を区別することができず、常に両方出力されます。 また、マッピング後にスタッツを表示してくれない点が、Bowtieと比べるとやや不便です。
colorspace read に関しては、 csfastqを作成すればマッピング可能と書いてあったりするのですが、自分は成功したことがありません。 うまくいくやり方をご存知の方はご教示ください。。

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

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

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

pyenvでPython環境を構築する

ざっくりですが。

Pythonツールのインストール

最近はPythonで書かれたツール類が増えてきています。 pip や conda でツールをインストールできるのはエラーも少なく、大変便利ですね。
ただ、pythonは2系と3系があり、python2.xでしか動かないツール、python3.xでしか動かないツールがあります。 違うバージョンで動かそうとするとエラーになったりして、ローカルPCなら適宜使い分ければいいのですが、 共有サーバ上で使い分けたりする場合はなかなかしんどいです。 （今そういう状態になっている）

ここでは、python2系と3系を使い分ける方法を簡単に記載します。
ここに書かれている内容は以下のエントリをベースにしていますので、 anaconda, pyenvとは何ぞや？という人はまず以下のエントリを参照ください。

データサイエンティストを目指す人のpython環境構築 2016

Win, MACのローカルPCの場合

Anacondaを直接インストールするのが一番簡単だと思います。
Anacondaも2系と3系がありますが、3系をインストールし、必要があれば中で2系の環境を構築するのがよいでしょう。
（構築方法は上記エントリを参照）

pyenvで環境構築

Linuxのターミナル上で動かす場合は、pyenvが一番良いのではないかと今のところ思っています。
システムにインストールされているpythonを色々いじるのはリスクが大きいことと、 pyenv内であれば、よくわからない状態になってしまった場合に新規入れ直しが簡単なためです。
anaconda3の中でpy2環境を構築するのでも良いのですが、一部のいけてないツールはpy2環境のbinを読みに行ってくれなかったりするので、pyenvで切り替える方がエラーが少なく済むかもしれません。
ただ、pyenvの環境設定は基本的にユーザ単位のようなので、共有サーバ上で使う場合、各ユーザが個別にpyenvをインストールし、設定作業をする必要があるようです。

pyenvのインストール

git がインストールされている前提です。

 $ git clone https://github.com/yyuu/pyenv.git ~/.pyenv
 $ echo 'export PYENV_ROOT="$HOME/.pyenv"' >> ~/.bashrc
 $ echo 'export PATH="$PYENV_ROOT/bin:$PATH"' >> ~/.bashrc
 $ echo 'eval "$(pyenv init -)"' >> ~/.bashrc
 $ source ~/.bashrc

pyenvでanacondaのインストール

ここではpyenv localを使う前提でanacondaの2系と3系をインストールしています。内部でpy2 環境を構築する場合はanaconda 3のみインストールします。

 $ pyenv install anaconda2-4.4.0
 $ pyenv install anaconda3-4.4.0
 $ pyenv rehash
 $ pyenv global anaconda3-4.4.0

ツール類のインストール

何しろnumpyとscipyは必ず必要になりますので最初に入れておきましょう。

 $ conda install -y numpy scipy

例えばDeepToolsをインストールしたい場合は以下のようにします。

 $ conda install -c bioconda deeptools

biocondaには色々とツールが入っているのですが、バージョンが古かったりエラーがあったりしますので、注意してください。

Python 2.x でしか動かないもの

• MACS2 及び 内部でMACS2を動かすツール（Mangoなど）
• Segway 2.0
• Hi-C pro
• RSeQC
• qiime

など

Python 3.x で動くもの

• DeepTools
• HiClib
• mageck-vispr
• Snakemake

など

ENA,DDBJからfastqを取得する（2019/9/6追記）

前回のエントリでは、SRAからfastqを取得する方法を紹介しました。

SRAからfastqを取得する (2019/07/16 追記あり) - Palmsonntagmorgen

一方、ENA (European Nucleotide Archive) やDDBJから直接fastqファイルをダウンロードすることも可能です。

ENA Browser
http://trace.ddbj.nig.ac.jp/dra/index_e.html

SRA, DDBJ, ENAのデータベースは互いに同期されており、基本的にどのサイトからでも同一のファイルを取得できます。
2017/7/13の時点では、ENAからは fastq.gz 形式、DDBJからは fastq.bz2 形式でダウンロード可能です。
SRAを使うと、sraのダウンロード -> fastq への変換 -> fastq.gz への圧縮と三段階必要であるのに比べ、 上記サイトを使えば fastq.gz のダウンロードのみになるので、簡便かつ高速です。 なお私的には.bz2形式よりも.gz形式の方が使いやすいので、ここではENAからのダウンロードを紹介します。
DDBJからのダウンロードについては下記のサイトに詳述されていますので、そちらを参考にしてください（適当）。

https://bi.biopapyrus.jp/rnaseq/fastq-data/download-fastq-ddbj.html

（2019/9/6追記：DDBJは2017年4月7日からSRAのミラーリングを停止しているそうです。SRAのデータのダウンロードにはSRA本体かENAを使ってください。
　参考：https://www.ddbj.nig.ac.jp/news/ja/180918_2.html）
（2021/5/7追記：2020年からDDBJのミラーリングは再開されたそうです。）

ENAからデータの取得

ここではサンプルデータとして SRR390728 を使います。
ENAのサイト上でSRR390728を検索すると、下記のような画面になります。

"FASTQ files (FTP)" のリンクをクリックするとダウンロードが始まります。 SRR390728 は paired-end なので、２ファイルに分かれています。 singleとpairedを間違うこともないので、簡単ですね。
Linux ターミナル上ではwgetを使ってダウンロードをすることも可能です。 上記ページからFTPのリンクをコピーし、以下のようにコマンドを打ちます。

 $ wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR390/SRR390728/SRR390728_1.fastq.gz
 $ wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR390/SRR390728/SRR390728_2.fastq.gz

wget を使ったファイルの一括取得

サンプル数が多い場合にひとつひとつ手作業でダウンロードするのは面倒です。 そのような時はFTPリンクをファイルに保存し、wget -i を使うと簡単です。

ここではサンプルデータとして SRP000712 を使います。
ENA上で SRP000712 を検索すると、以下のように25サンプル含まれていることがわかります。

画面上部の "TEXT" をクリックすると、各サンプルの詳細をテキスト形式でダウンロードできます。 この例では "PRJNA106023.txt" というファイル名で保存されます。
FTPのURLは10列目に記載されているので、10列目だけ抽出したファイルを作成しましょう。

 $ cut -f10 PRJNA106023.txt > ftp.txt

あとはこのファイルにwget -i を実行すれば順番にダウンロードされます。

 $ wget -i ftp.txt

サンプルによっては１サンプルにつき複数のURLが ; 区切りで記載されていたりしますが、 ; を改行に変換するなどして対応してください。

SRAからfastqを取得する (2019/07/16 追記あり)

更新の間が随分空いてしまいました。 その間に２つの学会に参加してきたのですが、海外の解析手法の進化具合にずいぶん衝撃を受けました。
が、ここは予定通り初歩的な作業から説明していきたいと思います。

今日はSRA(Sequence Read Archive) からfastqファイルを取得する方法です。

SRA Toolkit

まずはNCBIのサイトから SRA Toolkit をダウンロードします。

SRA Toolkit download

MAC用、Win用、Linux用と分かれているので、適切なものを選んでください。
例えばubuntu64なら、ダウンロードした圧縮ファイル(sratoolkit.2.8.2-1-ubuntu64.tar.gz)を解凍後、sratoolkit.2.8.2-1-ubuntu64/bin ディレクトリにパスを通せばOKです。
SRA Toolkitはかなり頻繁にアップデートされるので、ちょくちょく最新版をチェックしておくことをおすすめします。
2017/7/5時点での最新版は 2.8.2-1 です。

fastq-dumpでfastqをダウンロード (single-end)

ここではサンプルデータとして SRR390728 を使います。

 $ fastq-dump SRR390728

single-end データの場合、基本的にはこれだけでOKです。 Hiseqデータならば.fastq形式で、SOLiDでシーケンスされたcolor spaceデータならばcsfastq形式で出力されます。
圧縮形式(.fastq.gz)で出力したければ --gzip オプションを使います。

 $ fastq-dump SRR390728 --gzip

その他のオプションが知りたい場合は--help オプションを使用してください。

fastq-dumpでfastqをダウンロード (paired-end)

paired-endの場合、デフォルトで変換するとF3, F5リードがつながったSRR390728.fastq ファイルが出力されてしまいます。
F3, F5リードを分割して出力するには --split-files オプションを追加します。

 $ fastq-dump SRR390728 --split-files

これで、SRR390728_1.fastq, SRR390728_2.fastq と分割されたfastqファイルが生成されます。
なお、single-end データに --split-files オプションを使った場合、XXX_1.fastq ファイルのみ生成され、内容はオプションなしの場合と同一です。 ので、single-end とpaired-endデータが混在していて面倒という場合には、全てのファイルに --split-files オプションを付加しても問題ありません。

fastq-dumpで.sraから.fastqに変換

上記コマンドを実行した場合、fastq-dumpコマンドは内部で

(1) .sraファイルのダウンロード（~/ncbi/public/sraに保存される）
(2) .sraファイルの.fastqへの変換

を行っています。

fastq-dumpでダウンロードされた.sraファイルは.fastq変換完了後も削除されませんので、 何らかの理由で.fastqファイルへの変換に失敗した（paired-endなのに--split-files オプションをつけ忘れたなど）場合には、 ダウンロード済の.sraから.fastqの変換を行うと二度手間になりません。

 $ fastq-dump ~/ncbi/public/sra/SRR390728.sra

こうすると、 SRR390728.sra から SRR390728.fastq への変換のみが行われます。
逆に、.sraファイルをいつまでも保存していると容量を圧迫しますので、必要なくなった.sraファイルはまとめて削除しましょう。

(2019/07/16 追記) version 2.9.1 から、fasterq-dumpというコマンドが登場しています。fastq-dumpの後継で、マルチスレッドにすることで高速化しているようです。fastq-dumpは将来的に廃止されるそうですので、今後はこちらを使うとよいでしょう。なお--gzip オプションは消えているようですので、生成したfastqを pigzなどで圧縮するようにしましょう。

prefetchでsraをダウンロード

fastqファイルの取得・変換には時間がかかります。
大量のデータを一度に取得する場合には、 安全を期してまず.sraファイルを取得したいという場合もあるでしょう。
prefetchコマンドを使うと、sraファイルのダウンロードのみを行うことができます。

 $ prefetch SRR390728

このコマンドを実行すると、SRR390728.sra ファイルが ~/ncbi/public/sraに保存されます（出力先ディレクトリを何故指定できないのかは謎です…）。

--option-file コマンドを使うと、一度に複数の.sraファイルをダウンロードすることもできます。

 $ prefetch --option-file SRR_Acc_List.txt 

こうすると、SRR_Acc_List.txt に記載されているサンプル全てをダウンロードします。
これらのオプションはツールのバージョンアップと共に変わっていく可能性がありますので、新しいバージョンを使う場合には動作確認を忘れずに。

2bit genome を作成する

2bit genome はゲノム配列ファイルを2bit （バイナリ）形式で格納したものです。
2bit 形式はテキストエディタで開くことはできませんが、multifasta 形式よりも非常に高速にプログラムに読みこむことができるため、
ゲノム解析ツールを使う際にまれに 2bit genomeを要求されることがあります。
初めて見ると戸惑うかもしれませんが、UCSC genome browserで提供されているので、それをダウンロードして使えばOKです。
例えばhg38であれば http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/ の hg38.2bit がそれに該当します。

2bit genome の作成

この2bit genomeファイルはUCSCから提供されているプログラムで作成することも可能です。
UCSC genome browserからダウンロードした ゲノム配列データにはコンティグ配列も含まれていますので（下記エントリ参照）、
それらを取り除いた2bitファイルが欲しいというような場合は自分で作成するとよいでしょう。

rnakato.hatenablog.jp

twoBitToFa と faToTwoBit を使った相互変換

2bit 形式 と multifasta 形式の相互変換には twoBitToFa と faToTwoBit のふたつのプログラムを使います。
まずはこの２つのプログラムをダウンロードし、実行権限を付与します。

 $ wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
 $ wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/twoBitToFa
 $ chmod +x faToTwoBit twoBitToFa  # 実行権限の付与

UCSC genome browserからhg38 の 2bit genomeをダウンロードしてみましょう。

 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit
 $ less hg38.2bit   # lessコマンドでオープンを試みると
"hg38.2bit" may be a binary file.  See it anyway? 

hg38.2bitはバイナリファイルなので開けませんよ、と警告が出ます。

2bit genome → multifasta の変換

2bit genome → multifasta の変換には twoBitToFa を使います。

 $ ./twoBitToFa hg38.2bit hg38.fa
 $ less hg38.fa
>chr1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
...

multifastaはテキストファイルなので、lessでオープンすることができます。
なお、ここでNが並んでいるのは、chr1の末端はテロメア領域であり配列が解読されていないためです。 このような未解読の領域は「ギャップ領域」と呼ばれます。
ヒトゲノム解読が完了したとは言っても、このようなギャップ領域はまだたくさん残っているということですね。

multifasta → 2bit genome の変換

こちらは faToTwoBit を使います。

 $ ./faToTwoBit hg38.fa hg38.2bit

hg38.faとhg38.2bitはファイル形式が異なるだけで内容は同一ですので、
このように、何度でも相互変換をすることができます。

なお、http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/ に格納されている UCSCのコマンドツール群には他にも便利なものが色々ありますので、
全部ダウンロードして利用可能にしておくのがオススメです。