blasrのインストール（cmakeを使用する）

ずいぶん久しぶりの更新（100日振り！）。でも忘れているわけではないのです・・・
Linux（ubuntu 16.04 LTS)にblasrをインストールした。いろいろ苦戦したので、方法の簡単なメモを残しておく。
blasrについては、以前も別のコンピュータにインストールしことがあるが、その時といろいろ違っていたので、改めて記録しておく（またすぐに方法が変わる可能性もあるけれど・・・）。

今回は、前回とは異なり、cmakeを使う方法でインストールした。普通にmakeを使う方法は、自分の環境ではうまく行かなかった。

（１） 下準備
最初に、blasrはPython 2.7で動くので、LinuxのPythonが3系になっている場合は、2.7にしておく。以下を参照。
Phython3.6をインストールして切り替える - NGSデータ解析まとめ

cmakeをインストールする。以下のサイト参照にした。普通にapt-get installで問題なくインストールできた。
次に、Ninjaをインストールする。以下のコマンドによる。

sudo apt-get update
sudo apt-get install ninja-build

次に、BOOSTをインストール。このサイトを参照に、sudo apt-get install で普通にインストールできる。
さらに、htslibをインストールする。Githubから最新版をダウンロードしてコンパイルする。

（2） HDF5（libhdf5）のインストール
blasrをコンパイルするために、HDF5が必要なのでインストールする。また、Ubuntu 16.04の場合は、このサイトを参照に、libhdf5をインストールでいいと思う（たぶん）。

以下のようにする。

sudo apt-get install libhdf5-dev

（3） blasrのインストール

インストール方法は、以下のサイトを参考にした。
blasr/INSTALL_CMAKE.md at master · PacificBiosciences/blasr · GitHub

まず最初に、適当なpath（ここでは/home/hashi-linux2/）に"blasr_install"ディレクトリを作成し、gitでblasrをダウンロードする。ダウンロードした"blasr"ディレクトリに入る。

mkdir blasr_install
cd blasr_install
git clone https://github.com/PacificBiosciences/blasr.git --recursive
cd blasr

次に、gitでディレクトリ内容のアップデートを行う。

git submodule update --init --remote

ここで、cmake & ninjaによるcompile, buildを行う。

mkdir build && cd build    # build ディレクトリの作成、その中に移動
cmake -GNinja .. && ninja    # compile & build

うまくコンパイルが終了したら、ディレクトリ内に"blasr"のバイナリができる。ここでエラーが出る場合、メッセージにしたがって必要なライブラリをインストールしたりする必要がある。今回、特に引っかかったのは、このサイトにあるタイプのエラーで、このエラーが出た場合、"build"内にできたすべてのファイル、フォルダを削除して再度compileを行うとうまくいった。

ここまでできたら、今度は環境変数LD_LIBRARY_PATHを設定する。以下をターミナルに入力する。また、設定ファイル.profileにも書いておく（LD_LIBRARY_PATHは、それぞれの環境によって異なるので、適宜書き換えて下さい）。

export LD_LIBRARY_PATH="/home/hashi-linux2/blasr_install/blasr/libcpp/alignment:/home/hashi-linux2/blasr_install/blasr/libcpp/hdf:/home/hashi-linux2/blasr_install/blasr/libcpp/pbdata:/usr/lib/x86_64-linux-gnu/hdf5/serial/:/usr/local/lib/"
# メモ：最後の/usr/local/lib/は"libhts.so.2"の位置（環境によって場所は異なる）

動作確認のため、以下のコマンドを入力する。

./blasr --version

これでusageのメッセージが出れば、compileがうまくいっているので、以下のコマンドでインストールする。

sudo ninja install

これで終了！　やれやれ。

Python3.6をインストールして切り替える

NGS解析のとき、使用するLinux (Ubuntu 16.04)でPyhton 3.6とPyhton 2.7を切り替えることがあるので、その方法についてのメモ。
pyenvを使用する。

以下のサイトを参考にした。pyenvのインストールから、python3.6.0のインストールまではこちら。

Python3.6.0をUbuntu16.04に導入する。

使用するPythonのバージョンを3.6.0に切り替えるには、以下のようにする。

pyenv global 3.6.0

Pythonのバージョンを確認するには、以下。

python -V

ここでバージョンが切り替わっていれば、"Python3.6.0"と表示される。

再び、バージョンを2.7（システムのデフォルト設定）に戻すには、以下のようにする。

pyenv global system

使用できるPythonのリストは、以下で表示される。

pyenv versions

MacでフォーマットしたHDDをLinux (Ubuntu12.04)でマウントする（メモ）

HDDをフォーマットするときにしばしば迷うのでメモ。

MacでのHDDのフォーマット（MacOS拡張（ジャーナリング））は、Linux（Ubuntu）で普通に読み書き可能でマウントできるらしい。以下のURLによると、「MacOS拡張（ジャーナリング）」は、HFS+というフォーマット形式と同じとのこと。なるほど〜。

qiita.com

GATKを用いた多型サイトの抽出：-L optionについて

参考URL: https://www.broadinstitute.org/gatk/guide/article?id=4133

GATK HaplotypeCallerを使用して、リファレンスにマッピングされた特定の個人のヒトゲノムからSNPやINDELなどの多型を抽出するとき、-L optionを使用することで、計算量を節約できます。

たとえば、ヒトゲノムの20番染色体のみから多型を抽出する場合、-L optionの引数を以下のように指定します。

# 変数定義
sample=C8W82ANXX_PG1144_01A15_H1
echo ${sample}

# GATK HaplotypeCaller実行
java -Xmx4g -jar GenomeAnalysisTK.jar \
	-rf BadCigar -rf FailsVendorQualityCheck -rf MappingQualityUnavailable \
	-T HaplotypeCaller -R human_g1k_v37_decoy.fasta \
	-I ${sample}.fixmate.dedup.realign.recal.sort.addRG.bam \
	--dbsnp dbsnp_138.b37.vcf \
	--emitRefConfidence GVCF \
	--variant_index_type LINEAR \
	--variant_index_parameter 128000 \
        -L 20 \        # -L 調べたい染色体番号（scaffold, contig）を引数に指定
	-o ${sample}_raw_variants.20.g.vcf

また、たとえばexomeのみで多型を調べたい場合、Broad Instituteが提供している"Broad.human.exome.b37.interval_list"などを使用することで、exon領域のみの多型を抽出できます（"Broad.human.exome.b37.interval_list"についてはこちらを参照）。

# 変数定義
sample=C8W82ANXX_PG1144_01A15_H1
echo ${sample}

# GATK HaplotypeCaller実行
java -Xmx4g -jar GenomeAnalysisTK.jar \
	-rf BadCigar -rf FailsVendorQualityCheck -rf MappingQualityUnavailable \
	-T HaplotypeCaller -R human_g1k_v37_decoy.fasta \
	-I ${sample}.fixmate.dedup.realign.recal.sort.addRG.bam \
	--dbsnp dbsnp_138.b37.vcf \
	--emitRefConfidence GVCF \
	--variant_index_type LINEAR \
	--variant_index_parameter 128000 \
        -L Broad.human.exome.b37.interval_list \        # -L Broad.human.exome.b37.interval_listを引数に指定
	-o ${sample}_raw_variants.exome.g.vcf

もちろん、全ゲノム領域を対象にして多型を抽出したい場合は、-L optionなしで解析します。

BAM fileにRead Groupを付ける（GATKへの対応）

以下は個人的なメモ（覚え書き）になります。

GATKでは、BAM fileにRead Group (@RG)が付いてないとエラーが出て解析ができないようです。でもRead Groupって何、どうやって付けるの？　何を書けばいいの？　という点で疑問だったので、ちょっと勉強してみました。

まず、以下のページの日本語での解説によると、

http://blog.amelieff.jp/?eid=226752

BAM fileのheaderへのRead groupの入れ方としては、BWAの-R optionを使用するようです。またPicardを使って入れることもできます。実際には、両方やってみました。

## BWAの-R optionを使用する場合
# 変数の定義
ref=/PATH/to/human_g1k_v37_decoy.fasta
reads1=/PATH/to/QC.1.trimmed.fastq    # faQCs.plを使用した場合
reads2=/PATH/to/QC.2.trimmed.fastq
sample=C8W82ANXX_PG1144_02A16_H1_L003

#bwa //
bwa mem -t 12 -M \
-R "@RG\tID:FLOWCELLID\tSM:${sample}\tPL:illumina\tLB:${sample}_library_1" \
${ref} ${reads1} ${reads2} | samtools view -bS - > ${sample}.bam

ここで”C8W82ANXX_PG1144_02A16_H1_L003"というのは今回の自分のsample file (FASTQ)の名前で、特に何でもいいようです。とりあえずsample fileの名前にしています。

# picard AddOrReplaceReadGroupsを使用する場合

sample=C8W82ANXX_PG1144_02A16_H1_L003

java -jar picard.jar AddOrReplaceReadGroups \
    INPUT=${sample}.sort.bam \
    OUTPUT=${sample}.sort.addRG.bam \
    RGID=FLOWCELLID \
    RGLB= ${sample}_library1 \
    RGPU=H0164ALXX140820.2 \
    RGPL=illumina \
    RGSM=${sample}
# 181219追記：エラーが出るときは、以下のオプションをつける
    VALIDATION_STRINGENCY=LENIENT \
# 181219追記：BAM indexも作成
     CREATE_INDEX=true

Read group (@RG)をこのように付けることで、GATKはエラーを返さずに動くようになりました。

では、Read groupの各項目は何を意味しているのでしょうか？

Read groupの各項目の説明は、SAM/BAM fileについて説明したPDF（samtoolsのHPから）に書かれています。以下、GATKを使うときに入力が求められる5項目についての説明です。

項目	説明
ID	SAM/BAM fileの識別名（identifier, ID）。各@RG行は、固有のIDを付けなければならない。
LB	ライブラリ（何だろう？　ライブラリの通し番号？　上記のblogの例では"sample"と、別の本の例では"library1"など）
PU	Plathome unit（Illuminaの場合はflowcell-barcode-lane, Soildの場合はslide）
PL	配列を読むときに使用したプラットフォーム。有効（valid）な値はCAPILLARY, LS454, ILLUMINA, SOLID, HELICOS, IONTORRENT, ONT, PACBIOとのこと
SM	サンプルの名前

項目PL以外は、特にこの値でなければいけない、ということもなさそうなので、適当に入れれば良い？