2019 9/20追記 ゲノムのダウンロード
この記事では、初めてコマンドで動作するツールを使う方向けにゲノムの指定した領域から配列を抽出する方法について説明します。コンピュータはmacを想定しています。普通はpython3やanacondaを入れ、condaのコマンドを使ってBicoondaのプリビルドバイナリを導入するのが簡単ですが、$PATHの設定などでつまづくかもしれません。そこで、計算速度が遅くなることは気にせず、dokcerを使って仮想環境でBEDToolsを動かします。目的は、指定した領域の塩基配列を抽出し、エンリッチされた配列があるか調べることです。
１、macosのバージョンチェック
dockerをインストールするためにはmac10.12以上が必要です（間も無く10.13以上が最低ラインになるかもしれません）。また、メモリも8GBは必要です。できるだけ新しいmacを使ってください。
上のマシンはmacos10.12なのでインストール可能。メモリは96GBになっている。ストレージ（HDDやSSD）とメモリを混同しないように注意。
２、dockerのインストール
Qiitaの記事を参考に、dockerをインストールします。
dockerインストールに成功したら、右上にクジラのアイコンが表示されているはずです。
クジラのアイコンをクリックし、dockerの仮想環境で利用するメモリやCPU数を指定します。
mac proはCPUやメモリのリソースが多いため、CPU18、メモリ84GBまで増やしています。 下のボタンが緑色になったら使用できます。
３、ターミナルを起動します。アプリケーション（commnad + shift + A）=>ユーティリティにあります。
起動したら、dockerと打ちます。
docker
ヘルプが出ればインストールに成功しています。
４A 、ヒトゲノムアセンブリのダウンロード
Homo sapiens - Ensembl genome browser 97
右の方のDownload FASTA を選択。ウィンドウが出てくるのでゲストを選択。
dnaを選択。全ゲノムならばHomo_sapiens.GRCh38.dna.toplevel.fa.gzをダウンロードする。rmとsmはリピートがハードマスクかソフトマスクされたアセンブリになる。
解凍しておく。
またはNCBIその他からもダウンロードできる。
https://www.ncbi.nlm.nih.gov/genome/51
genomeをクリックする。
４B 、データを作業フォルダに集める。
ここではデスクトップにhumanというフォルダを作り、そこにヒトのGRCh38アセンブリをコピーしました。
ターミナルで移動します。ターミナルでcd （change directory）と打ち、それからスペースキーを一度打って１スペース空け、humanのフォルダをターミナルにドラッグ&ドロップします。
こんな感じになるので、そのままエンターキーを押したら移動します。
ls（list segments）+Enterと打って用意したファイルがあるか確認します。
ls
５、bedtoolsを含むイメージのpull
dokcerhubというイメージ共有サイトからbedtoolsがプリインストールされたイメージをローカルにダウンロードします。ここではlindsayliangさんのイメージを使います。以下のコマンドをターミナルにコピペしてEnterして下さい。
docker pull lindsayliang/bedtools
６、イメージの動作確認
全コンテナのpullが終わったら helpが表示されるか確認します。以下のコマンドをコピペして下さい。
docker run --rm -itv $(PWD):/data -w /data lindsayliang/bedtools bedtools
これはローカルと仮想環境の共有ディレクトリをカレントディレクトリ（現在いるフォルダ）、run後の初期パスをカレントディレクトリ、exit後に消す設定でイメージをrunし、イメージにインストールされたbedtoolsのヘルプを表示しています。
うまく動けば以下のようなbedtoolsのhelpメッセージが表示されるはずです。
$ docker run --rm -itv $(PWD):/data -w /data lindsayliang/bedtools bedtools
bedtools: flexible tools for genome arithmetic and DNA sequence analysis.
usage: bedtools <subcommand> [options]
The bedtools sub-commands include:
[ Genome arithmetic ]
intersect Find overlapping intervals in various ways.
window Find overlapping intervals within a window around an interval.
closest Find the closest, potentially non-overlapping interval.
coverage Compute the coverage over defined intervals.
map Apply a function to a column for each overlapping interval.
genomecov Compute the coverage over an entire genome.
merge Combine overlapping/nearby intervals into a single interval.
cluster Cluster (but don't merge) overlapping/nearby intervals.
complement Extract intervals _not_ represented by an interval file.
shift Adjust the position of intervals.
subtract Remove intervals based on overlaps b/w two files.
slop Adjust the size of intervals.
flank Create new intervals from the flanks of existing intervals.
sort Order the intervals in a file.
random Generate random intervals in a genome.
shuffle Randomly redistrubute intervals in a genome.
sample Sample random records from file using reservoir sampling.
spacing Report the gap lengths between intervals in a file.
annotate Annotate coverage of features from multiple files.
（以下略）
７、指定した領域のBEDファイルの準備
ここからが本番です。取ってきたい領域を指定したファイルを用意します。chr1の三箇所、chr3の一箇所から取ってくるには以下のようにします。excelを使っています。
左端に染色体番号、２列目にスタート（上流側）ポジション、３列目にエンド（下流側）ポジションを記入。
Miなどのテキストエディタに配列をコピペ。region.bedというファイル名で、作成したフォルダに保存する（commnad + S）。
右上の改行コードボタンの表示がCR（キャリッジリターン）になっていたら、LF（ラインフィード）に変更。
８、配列の抽出
作業前にindexファイルを作成します。これは本の索引に相当するファイルで、プログラムが高速動作するためにゲノムのファイルと同じ場所にある必要があります。
docker run --rm -itv $(PWD):/data -w /data lindsayliang/bedtools \
samtools faidx Homo_sapiens.GRCh38.dna.toplevel.fa
ゲノムが大きいので10分以上かかります。 取ってきたい領域が特定の染色体に限定されるなら"-r xxx"という風にindex領域を指定することで時間短縮できますが、ここでは全領域のindexファイルを作成します（コマンドに慣れているならrefgenieというツールも使えます）。 Homo_sapiens.GRCh38.dna.toplevel.fa.faiが出力されます。
bedtools（紹介）を動かします。途中にある\はコマンドの途中に改行する時に使います。
docker run --rm -itv $(PWD):/data -w /data lindsayliang/bedtools \
bedtools getfasta -fi Homo_sapiens.GRCh38.dna.toplevel.fa \
-bed region.bed > reads.fasta
ファイル名が違うところは変えて下さい。
reads.fastaが出力されます。lessコマンドで中身を見てみます。
less reads.fasta
出力されてますね。スペースキーで１ページ下にスクロールします。bで１ページ戻ります。qを押すとlessコマンドは終了します。
exitでターミナルは終了しますが、commnad + Qでターミナルをアプリごと閉じてもらっても構いません。
追記
取り出す領域が１つならsamtools faidxが早いです。
docker run --rm -itv $(PWD):/data -w /data lindsayliang/bedtools \
samtools faidx Homo_sapiens.GRCh38.dna.toplevel.fa chr1:10000-10010 \
> reads.fasta
９、モチーフ検索
MEME CHIPにアクセスして、エンリッチされているモチーフがあるか調べます。
http://meme-suite.org/tools/meme-chip
先ほど出力された.fastaファイルをアップロードします。
MEME CHIPについてはNGS_handson2015の説明を参考にして下さい。
https://github.com/suimye/NGS_handson2015/wiki/PeakCallAndMDA
関連
