macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

初めてコマンドを使う人向けの解説:その1、指定した領域から配列を抽出する

2019 9/20追記 ゲノムのダウンロード

 

この記事では、初めてコマンドで動作するツールを使う方向けにゲノムの指定した領域から配列を抽出する方法について説明します。コンピュータはmacを想定しています。普通はpython3やanacondaを入れ、condaのコマンドを使ってBicoondaのプリビルドバイナリを導入するのが簡単ですが、$PATHの設定などでつまづくかもしれません。そこで、計算速度が遅くなることは気にせず、dokcerを使って仮想環境でBEDToolsを動かします。目的は、指定した領域の塩基配列を抽出し、エンリッチされた配列があるか調べることです。

 

1、macosのバージョンチェック

dockerをインストールするためにはmac10.12以上が必要です(間も無く10.13以上が最低ラインになるかもしれません)。また、メモリも8GBは必要です。できるだけ新しいmacを使ってください。

f:id:kazumaxneo:20190919161537j:plain

上のマシンはmacos10.12なのでインストール可能。メモリは96GBになっている。ストレージ(HDDやSSD)とメモリを混同しないように注意。

 

 

、dockerのインストール

Qiitaの記事を参考に、dockerをインストールします。


dockerインストールに成功したら、右上にクジラのアイコンが表示されているはずです。

f:id:kazumaxneo:20190919163339j:plain


クジラのアイコンをクリックし、dockerの仮想環境で利用するメモリやCPU数を指定します。

f:id:kazumaxneo:20190919162213j:plain

mac proはCPUやメモリのリソースが多いため、CPU18、メモリ84GBまで増やしています。 下のボタンが緑色になったら使用できます。

 

 

、ターミナルを起動します。アプリケーション(commnad + shift + A)=>ユーティリティにあります。

f:id:kazumaxneo:20190919151618p:plain

起動したら、dockerと打ちます。

docker

 ヘルプが出ればインストールに成功しています。

f:id:kazumaxneo:20190919151713p:plain

 

4A 、ヒトゲノムアセンブリのダウンロード

Homo sapiens - Ensembl genome browser 97

f:id:kazumaxneo:20190920133334p:plain

右の方のDownload FASTA を選択。ウィンドウが出てくるのでゲストを選択。

f:id:kazumaxneo:20190920133424p:plain

dnaを選択。全ゲノムならばHomo_sapiens.GRCh38.dna.toplevel.fa.gzをダウンロードする。rmとsmはリピートがハードマスクかソフトマスクされたアセンブリになる。

f:id:kazumaxneo:20190920133713p:plain

解凍しておく。
 

またはNCBIその他からもダウンロードできる。

NCBI

https://www.ncbi.nlm.nih.gov/genome/51

f:id:kazumaxneo:20190920135434p:plain

genomeをクリックする。

 

 

4B 、データを作業フォルダに集める。

ここではデスクトップにhumanというフォルダを作り、そこにヒトのGRCh38アセンブリをコピーしました。

f:id:kazumaxneo:20190919164145j:plain

 

ターミナルで移動します。ターミナルでcd (change directory)と打ち、それからスペースキーを一度打って1スペース空け、humanのフォルダをターミナルにドラッグ&ドロップします。

f:id:kazumaxneo:20190919164803j:plain

 

こんな感じになるので、そのままエンターキーを押したら移動します。

f:id:kazumaxneo:20190919164940j:plain

 

ls(list segments)+Enterと打って用意したファイルがあるか確認します。

ls

 

 

、bedtoolsを含むイメージのpull

 dokcerhubというイメージ共有サイトからbedtoolsがプリインストールされたイメージをローカルにダウンロードします。ここではlindsayliangさんのイメージを使います。以下のコマンドをターミナルにコピペしてEnterして下さい。

docker pull lindsayliang/bedtools

f:id:kazumaxneo:20190919153942p:plain

 

 

、イメージの動作確認

全コンテナの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を使っています。

f:id:kazumaxneo:20190919183006j:plain

 左端に染色体番号、2列目にスタート(上流側)ポジション、3列目にエンド(下流側)ポジションを記入。

 

Miなどのテキストエディタに配列をコピペ。region.bedというファイル名で、作成したフォルダに保存する(commnad + S)。

f:id:kazumaxneo:20190919183425j:plain

右上の改行コードボタンの表示が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

f:id:kazumaxneo:20190919185437j:plain

出力されてますね。スペースキーで1ページ下にスクロールします。bで1ページ戻ります。qを押すとlessコマンドは終了します。 

 exitでターミナルは終了しますが、commnad + Qでターミナルをアプリごと閉じてもらっても構いません。 

 追記

取り出す領域が1つなら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

f:id:kazumaxneo:20190919201134j:plain

 

先ほど出力された.fastaファイルをアップロードします。

f:id:kazumaxneo:20190919201937j:plain

 

MEME CHIPについてはNGS_handson2015の説明を参考にして下さい。

https://github.com/suimye/NGS_handson2015/wiki/PeakCallAndMDA

 

関連