ゼロからはじめるChIP-seq解析

# はじめに

ChIP-seq(chromatin immunoprecipitation followed by sequencing)は特定の転写因子の結合やヒストン修飾がゲノム上のどの位置でどれぐらいの頻度で起こっているのかを網羅的に測定する方法です. この記事を読めば環境構築から始めて論文のデータを用いて解析を行い, 最終的に結果をゲノムブラウザで見たり, ピークコールをやったりするところまでできるようになりますよ～. それでは早速環境構築から始めましょう！

# 環境構築

Windows を使っている方は Windows Subsystem for Linux (WSL) をインストールして以下の Linux の手順を実行してください. また, WSL 最大のハマりどころについて下のヘルプ4で解説しています.  

## miniconda3のインストール

まずはパッケージマネージャーであるminicondaをインストールします. パッケージマネージャーは環境構築時のツールのインストールや管理を簡単にできるようにしてくれます.  
  
Linuxの場合

```
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

bash Miniconda3-latest-Linux-x86_64.sh
```

Macの場合

```
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh

bash Miniconda3-latest-MacOSX-x86_64.sh
```

これを実行したら指示に従ってENTERを押したり, `yes`と入力したりしてください. 質問にはすべて`yes`と入力していただいて構いません. すべて終わったら一度ターミナルを閉じてください. また起動させて次はminicondaにチャンネルを追加するので以下を実行してください. __必ずこの順で実行してください__.

```
conda config --add channels defaults
conda config --add channels conda-forge
conda config --add channels bioconda
```

ここでcondaがないよと言われてしまった人は下のヘルプ1を参照してください.  

## 必要なツールのインストール

それではminicondaを使って必要なツールをインストールしていきます. 今回は以下のツールをインストールします.  

- [SRA Toolkit](https://github.com/ncbi/sra-tools) SRAからシークエンスデータをダウンロードするのにつかうツール
- [Trimmomatic](http://www.usadellab.org/cms/?page=trimmomatic) アダプタートリミングをするツール
- [fastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) fastqファイルのQCを行うツール
- [Bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) マッピングをするツール
- [Picard](http://broadinstitute.github.io/picard/) PCR duplicateを除去するツール
- [samtools](http://samtools.sourceforge.net/) 色々なところで使うツール 
- [deepTools](https://deeptools.readthedocs.io/en/develop/) データの補正に使うツール
- [HOMER](http://homer.ucsd.edu/homer/ngs/) 色々できるツール


以下のようにminicondaの`conda`コマンドを使ってインストールします.  

```
conda install sra-tools
conda install trimmomatic
conda install fastqc
conda install bowtie2
conda install picard
conda install samtools
conda install deeptools
conda install homer
```

`Proceed ([y]/n)?`にはすべて`y`と答えてください. 
  
これで環境は整いました！これからデータを触っていきますが本記事では混乱を避けるために以下**すべてのコマンドを同一のディレクトリで行った場合を想定して書いています**. それではデータをもらいに行きましょう～.

# シークエンスデータのダウンロード

## データのありかを確認する

[SRA(Sequence Read Archive)](https://www.ncbi.nlm.nih.gov/sra) からシークエンスデータをダウンロードします. 今回は[Kagey MH *et al.,*](https://www.nature.com/articles/nature09380)(GSM560348)からマウスのES細胞のMed1のChIP-seqのデータと[Creyghton *et al.,*](https://www.pnas.org/content/107/50/21931.long)(GSM594579)から同じくマウスのES細胞のH3K27AcのChIP-seqのデータを使用します. また, Med1と同じ論文からインプットのシークエンスデータ(GSM560357)も取得します. ()内の番号は各データのGEO Accession番号です. NCBIのGene Expression Omnibus(GEO)でこれらのデータのページ飛ぶのに必要な番号です.  

さて, データをダウンロードする前にまずはこれらのデータのありかをブラウザ上で確認しておきましょう. Med1のデータを例に説明します. まずは[GEOのサイト](https://www.ncbi.nlm.nih.gov/geo/)をブラウザで開いてください. 下の画像で赤く囲っている検索窓に`GSM560348`と打ち込んでください.

![GEO1.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/0897f373-0e7c-5c96-75b5-79d412396413.png)
  

下の画像のようにこのデータについての色々な情報が書いてあるのでしっかり目をとおしてください.

![GEO2.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/40bf9c1e-1872-a869-417a-446db971722c.png)


では, このページの下の方にあるSRAと書かれたところの右にある番号をクリックしてください.

![GEO3.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/427e31df-8e2a-c437-082f-c2962d04e763.png)


すると以下のようなページに行きます. 画像で赤く囲っている番号がSRR番号と言ってこのデータのアクセス番号でダウンロードする際に必要なものです.

![GEO4.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/ba16dca8-8471-e787-3de0-2445e219e78b.png)


論文中にはこれらの番号がどこかに書いてあってデータのありかがわかるようになっています.  

## SRAからデータをダウンロードする

sratoolkitの`fastq-dump`というコマンドを使用します. 使い方はとっても簡単でデータがシングルエンドの場合は

```
fastq-dump ダウンロードしたいデータのSRR番号 —-gzip
```

とするだけです. `—-gzip`オプションはデータをgzipに圧縮して落とすためのオプションです. 通信環境によってはこれを付けないとデータが重すぎてよく失敗するのでつけることをお勧めします.  

ペアエンドの場合は

```
fastq-dump ダウンロードしたいデータのSRR番号 —-gzip --split-files
```

としてください.   

これでこのコマンドを実行したディレクトリにシークエンサーから出力される生データであるfastqファイルがダウンロードされるはずです.  

それでは今回のデータをダウンロードします.

```
fastq-dump SRR058988 —-gzip #Med1
fastq-dump SRR066767 —-gzip #H3K27Ac
fastq-dump SRR058997 —-gzip #インプット
```

以下のように書き並べて一回でダウンロードすることもできます.

```
fastq-dump SRR058988 SRR066767 SRR058997 —-gzip
```

この処理ではSRA特有の圧縮ファイル形式である.sraファイルをダウンロードしてそれをfastqファイルに変換しています. 時間がかかるので気長に待ちましょう. 


ダウンロードが終わったらまず`gunzip`コマンドでファイルを展開します.  

```
gunzip SRR058988.fastq.gz
gunzip SRR066767.fastq.gz
gunzip SRR058997.fastq.gz
```  

が, `fastq-dump` は遅すぎるし挙句の果てには失敗したりするので NCBI が新たに開発した `fasterq-dump` を使うのがおすすめです.  

```  
fasterq-dump SRR058988 -e 8 -p #Med1
fasterq-dump SRR066767 -e 8 -p #H3K27Ac
fasterq-dump SRR058997 -e 8 -p #インプット  
```  

`-e` オプションで使用するコア数を指定します. また, `-p` をつけることで進捗を表示させることができるのであとどれくらいで終わるのか見えて心理的に楽です. なお, `fasterq-dump` はペアエンドかどうか自動で認識するので `--split-files` の指定は必要ありません.  


ダウンロードできたら, ファイル名がSRR番号で付けられているので名前を付け替えてわかりやすくしておきます.

```
mv SRR058988.fastq Med1.fastq
mv SRR066767.fastq H3K27Ac.fastq
mv SRR058997.fastq Input.fastq
```

# アダプタートリミング

[2021/3/7追記]
どうやら最近のmapperはトリミングをしなくてもよきに計らってくれるようなのでアダプターのトリミングはいらないかもです．トリムするアダプターは指定せず，クオリティーが低いリードを取り除くという目的でtrimmomaticを使うのはありだと思います．  

（参考）[When should I trim my Illumina reads and how should I do it?](https://dnatech.genomecenter.ucdavis.edu/faqs/when-should-i-trim-my-illumina-reads-and-how-should-i-do-it/)
  
  

それではシークエンス結果のクリーニングをTrimmomaticを使って行います. 

```
trimmomatic SE -threads 4 -phred33 Med1.fastq Med1_trimmed.fastq ILLUMINACLIP:$HOME/miniconda3/share/trimmomatic/adapters/TruSeq3-SE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
```

長いですが改行による事故を防ぐためにあえて改行していません. `trimmomatic`の直後には解析するデータがシングルエンド(SE)かペアエンド(PE)かを指定します. 次の`-threads`で使用するスレッドの数を指定します. `-phred33`はおまじないです. 必ず入力してください. その後にトリミングするファイルの名前, トリミング後のファイルの名前を続けて入力します.  

ILLUMINACLIPの後ろに書いてあるのはアダプター配列の情報の場所でminiconda3ディレクトリ以下にこのような場所にあるはずです. みなさんのminiconda3の場所に合わせて書き換えてください.（インストールする時に特に何もしていなければこのようにホームディレクトリにあるはずです.） また2:30:10はそれぞれ許容ミスマッチ数, palindrome clip threshold, simple clip thresholdを表します. 基本的にはここはいじらなくても大丈夫だと思います. また, `LEADING:3`, `TRAILING:3`はそれぞれリードの先頭, 末尾からクオリティスコアが3未満の塩基を削除するという意味です. `SLIDINGWINDOW:4:15`というのは4 bpごとにクオリティスコアの平均値を見ていき, それが15未満である部分を削除するという意味です. そして最後の`MINLEN:36`はリードの長さが36未満であるものを解析から外すという意味です. これらの設定は[Trimmomaticのページ](http://www.usadellab.org/cms/?page=trimmomatic)の"Quick start"のものを使用しました.
終わるとMed1_trimmed.fastqというファイルが生成されます. 残りの2つデータについても同じオプションで実行してください.

# トリミング後のfastqファイルのQC

トリミング後のfastqファイルのクオリティコントロールをfastQCを使って行います.  

```
fastqc --nogroup -o . Med1_trimmed.fastq
```

`--nogroup`を書いておくと3'末端のリードについても解析されます. `-o`に結果を出力するディレクトリを書きます. 最後にQCするファイルの名前を書きます.  

`-o`に指定したディレクトリに`Med1_trimmed_fastqc.html`というファイルができるのでブラウザで開いてみましょう.  

![fastqc.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/c354cfad-2864-b917-02fa-7ce8d6c7bdef.png)


左のSummaryのところがほとんど緑のチェックになっていたらクオリティはOKです. このデータきれいすぎる... 各指標については https://bi.biopapyrus.jp/rnaseq/qc/fastqc.html で説明してくださっています. こちらのサイト, bioinformaticsさんには常日頃からお世話になっています. いつもありがとうございます.  

残りの2つのデータについても同じオプションで行ってくださいね.  

# マッピング

いよいよBowtie2を使ってマッピングして行きます！と, その前にゲノムのインデックスをビルドしなくてはなりません. つまり, マッピングに必要なリファレンスゲノムの準備を行う必要があります.  

[UCSCのページ](https://hgdownload.soe.ucsc.edu/downloads.html)からmm10の全配列を`wget`コマンドでダウンロードします. `ref_genome`というフォルダを作ってそこに落とすことにします.

```
mkdir ref_genome
cd ref_genome
wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz
tar -xzf chromFa.tar.gz
```

また, 今回はrandomやunknownの配列やミトコンドリアゲノム(chrM.fa)は捨てることにします.

```
rm *random.fa
rm chrUn*
rm chrM.fa
```

残ったものを`cat`で一つの`mm10.fa`というファイルにします.

```
cat *.fa > mm10.fa
```

`ref_genome`と同じ階層にmm10_indexというディレクトリを作ってそこにインデックスをセーブすることにします.

```
cd .. #今ref_genomeにいる場合
mkdir mm10_index
bowtie2-build --threads 4 -f ./ref_genome/mm10.fa ./mm10_index/mm10_index
```

`bowtie2-build`がBowtie2のインデックスを作るコマンドです. `-f` オプションにゲノム配列の場所のパスを書いてその次にインデックスへのパスを書きます. `--threads`でスレッド数を指定できます. この場合, 先程作った`mm10_index`というディレクトリに`mm10_index*.bt2`というファイルが6つ作られるはずです.この操作は一回だけ行えばよいです. 

とても時間がかかります. この時間を使ってバイオインフォマティクスでの統計解析によく使われているR言語の勉強をしておきましょう.（いいねよろしくお願いします！）

https://qiita.com/roadricefield/items/001c882f84dd093f4407  

この記事ではRは使いませんが...  


※本当に時間がかかります. 帰る直前にコマンドを叩いて明日来たら続きをやるぐらいがいいレベルです. 僕も今日はここまでにしてアニメ観て寝ます.  

.........................  

みなさん, おはようございます！僕は`--threads`を指定するのを忘れてたので7時間もかかりました！それではマッピングします.

```
bowtie2 -p 4 -x ./mm10_index/mm10_index -U Med1_trimmed.fastq -S Med1.sam
```

`-p`はスレッド数, `-x`はインデックス, `-U`はマッピングするfastqファイル, `-S`は出力ファイル名です. 一旦samファイルとして出力します. ちなみにマッピングもかなり時間がかかります. 30分ぐらいですかね.  

終わったらsamtoolsでsamファイルをbamファイルに変換します.  

```
samtools view -b -o Med1.bam Med1.sam
```

これでマッピングは終了です！残りの2つデータについても同じオプションで行ってくださいね.

# PCR duplicate の除去（オプション）

最後にPicardを使ってPCR duplicateの除去を行います. 必ずしもやる必要はありません. 

```
samtools sort Med1.bam > Med1_sorted.bam #Picardを使うためにはbamファイルをソートする必要があります.

picard MarkDuplicates I=Med1_sorted.bam O=Med1_rm_dups.bam M=Med1_report.txt REMOVE_DUPLICATES=true
```

`I`にPCR duplicateを除去したいbamファイルの名前, `O`にduplicateを除去した後のファイルの名前, 計算結果のまとめであるレポートを作ってくれるのでにそれの名前を`M`に書きます.  

他の2つのデータも同じオプションで行ってください.  

そろそろ作業しているディレクトリがごった返してきたところなのでここでデータを整理します. Med1 ChIP-seqのデータは`Med1_data`というディレクトリに, H3K27Ac ChIP-seqのデータは`H3K27Ac_data`というディレクトリに, インプットのデータは`Input_data`というディレクトリに移しましょう.

```
mkdir Med1_data
mv Med1* Med1_data

mkdir H3K27Ac_data
mv H3K27Ac* H3K27Ac_data

mkdir Input_data
mv Input* Input_data
```

# ゲノムブラウザで結果を観察する

## bigWig ファイルを作成する

まずはbamファイルをゲノムブラウザで見やすいbigWigという形式に変換しましょう. これにはdeepToolsの`bamCoverage`を使用します. bigWigフォーマットでは決められたbin幅でのリード数の補正値が全ゲノムで計算されます. つまり, ゲノム上の各binでのChIP-seqのsignal強度が計算されます. これを作るためにはまずbamのインデックスファイルであるbam.baiファイルを作成する必要があるのでそれをsamtoolsで作りましょう.

```
samtools index Med1_data/Med1_rm_dups.bam
```

これを実行すると`Med1_rm_dups.bam.bai`というインデックスファイルが作成されます. 

では`bamCoverage`を実行します. 必ずbamファイルとそれのbam.baiファイルを同じディレクトリに置いて実行してください.

```
bamCoverage -b Med1_data/Med1_rm_dups.bam -p 4 --normalizeUsing RPGC --effectiveGenomeSize 2652783500 --binSize 1 -o Med1_data/Med1.bigwig
```
`-b`にbigWigファイルに変換するbamファイルの名前を書きます. `-p`はスレッド数. `--normalizeUsing`に各binで計算する補正値の種類を選びます. `RPKM`, `CPM`, `BPM`, `RPGC`, `None`が選べます. `None`を選んだ場合はbinに含まれるリード数がそのbinでの値になります. `--effectiveGenomeSize`にはゲノムのうちマッピング可能な部分の長さ(bp)を入力します. mm10(別名 GRCm38)の場合は`2652783500`です. (参考 https://deeptools.readthedocs.io/en/latest/content/feature/effectiveGenomeSize.html) `--binSize`には計算に使うbinの長さ(bp)を入力します. `-o`には出力ファイルの名前を書きます.

計算には時間がかかるのでその間にゲノムブラウザをインストールしましょう.

## IGV (Integrative Genomics Viewer)のインストール

ゲノムブラウザとはシークエンス結果を可視化するツールです. よく発表でゲノム上の特定の位置での〇〇-seqのシグナルが可視化されているのを見ますよね？あれです. それでは早速インストールしましょう！  

[IGVのダウンロードページ](https://software.broadinstitute.org/software/igv/download)で自分のOSのインストーラーをダウンロードしてください. IGVはGUI（グラフィックが出てきてマウスとキーボードで操作するインターフェイス）ですのでWSLを使っているWindowsユーザーの方はここではWindows版を選んでください. ダウンロードしたインストーラーを起動して指示に従ってインストールしてください. するとデスクトップに以下のようなIGVのショートカットが作られます. 

![IGV1.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/bd3c8a57-c96e-f1c9-255e-165de62c9d4c.png)
  

これをダブルクリックして起動させてください(起動に30秒ほどかかります). 起動後以下のようなウインドウが現れます.  

![IGV2.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/558efcf8-cfc4-96ac-c6b7-67b7e482b544.png)


今, hg19がロードされているのでmm10をダウンロードしてロードしましょう. 上の画面で赤く囲っているところの下矢印をクリックすると"More..."と出てくるのでそれをクリックしてください. そして"Mouse mm10"をクリックして左下の"Download Sequence"にチェックをいれて"OK"をクリック. これでmm10のダウンロードが始まります. 

![IGV3.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/138e1746-c4a6-9f2f-abbd-90cf22d53fd0.png)



そろそろ`bamCoverage`が終わったころですかね...？  

終わったら`Med1.bigwig`をIGVのウィンドウにドラッグアンドドロップしてください. 

![IGV4.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/7e3de502-d572-1050-b56e-af8a26a5d76d.png)


上の画像のようにMed1 ChIP-seqのプロファイルが表示されたでしょうか? この状態ではすべての染色体を俯瞰して見ており, 詳細がわからないので赤く囲ってある検索窓に色々な遺伝子名前を入力してそのgene bodyの場所に飛んでみましょう. 以下に一つだけ例を示します.

![IGV_Klf4.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/e18a9200-49d3-05cf-c7a5-65f31c0232e0.png)

*Klf4*は幹細胞性維持に重要な遺伝子ですがこのようにマウスのES細胞では転写を活性化する蛋白であるMed1が結合していることがわかります. こちらの画像では見やすいようにRefseq genesのトラックを引き上げています. また, 表示の縦軸のスケールがおかしい場合は"Med1.bigwig"という名前が書いてあるところで右クリックして"Autoscale"をクリックしてください.  

残りの2つのデータについても同様にbigwigを作成してIGVで確認してみてください.  

# ピークコール

それではシグナルのピークを統計学的基準に基づいて検出するピークコールを行います. 今回はHOMERの`findPeaks`を使います. これ以外にピークコーラーとしてよく使われているものに[MACS2](https://github.com/taoliu/MACS)があります. 興味があれば結果を比べてみてください.  

さて`findPeaks`を行う前にbamファイルをHOMERが扱うことのできるTagDirectoryという形に変換しなければなりません. これにはHOMERの`makeTagDirectory`を使います.  

```
makeTagDirectory Med1_data/Med1_TagDir -single Med1_data/Med1_rm_dups.bam
```

`makeTagDirectory`の直後に作成するTagDirectoryの名前を書き, その後にオプションを書き並べて最後にbamファイルの名前を書きます. 今回オプションはTagDirectory内をすっきりさせる`-single`オプションのみを指定しました. その他のオプションについては http://homer.ucsd.edu/homer/ngs/tagDir.html を参照してください. 残りの2つのデータについても同じようにTagDirectoryを作成してください.  

それでは`findPeaks`を走らせます.  

```
findPeaks Med1_data/Med1_TagDir -style factor -o auto -i Input_data/Input_TagDir
```

`findPeaks`の直後にピークコールを行うTagDirectoryの名前を書きます. `-style`オプションには今回はMed1を転写因子(transcription factor)とみなして`factor`と入力しました. `-o`オプションには結果の書き出し場所を書きますが`auto`としておけばピークコールを行うTagDirectory内に保存されます. また, `-i`オプションにはインプットのデータのTagDirectoryの名前を書きます. インプットのデータが無い時は`-i`オプション自体を入力しなければインプットなしで計算を行ってくれます.  

つぎにH3K27Ac ChIP-seqのデータのピークコールを行います.  

```
findPeaks H3K27Ac_data/H3K27Ac_TagDir -style histone -o auto -i Input_data/Input_TagDir
```

H3K27Acはヒストン修飾のデータなので`-style`オプションは`histone`としてください.  

それではピークコールの結果のファイルを見てみましょう. TagDirectoryの中に`peaks.txt`というファイルができていると思います(H3K27Acのほうは`regions.txt`という名前です). タブ区切りのテキストなのでExcelで開くと見やすいです.  


![findPeaks.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/044d2933-eaa7-8e5b-2e72-e8fcf86d9cbe.png)


上の方に様々な情報が書かれている後で画像の下半分のようなピークの情報が書かれています. `chr`, `start`, `end`の列を見れば各ピークのゲノム上での位置がわかります. `Normalized Tag Count`が各ピークのシグナル強度です. 10で割ればrpm(Reads per Million)になります. `findPeaks`の機能についてさらに詳しくは http://homer.ucsd.edu/homer/ngs/peaks.html を御覧ください.

次にこの結果を先程のbigwigファイルと合わせてIGVで見てみましょう. IGVにピークの情報を簡単に読み込ませるにはbedファイルという形式にするのがよいです. bedファイルは至ってシンプルで一番左から`chr`, `start`, `end`の3行で構成されているピークの位置情報を記録した以下の画像のようなテキストファイルです. (4行目以降にさらに情報が入っている場合もあります. ) 

![bed.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/cee9e82b-8125-08e7-cb16-cf749c502395.png)


`findPeaks`の出力の`peaks.txt`からは以下のコマンドで作成できます.

```
sed '/^#/d' Med1_data/Med1_TagDir/peaks.txt | awk -v 'OFS=\t' '{print $2, $3, $4}' > Med1_data/Med1_TagDir/peaks.bed
```

今作った`peaks.bed`を先程のIGVにドラッグアンドドロップしてみてください.  

![peak.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/6e48bd58-d569-092e-73a0-5399eec51cce.png)


このようにピークコールされた位置が表示されるはずです. わかりやすい...

# 結合モチーフ解析

転写因子はある程度決まった配列（結合モチーフ配列）に結合するのでピーク領域が与えられたとき, それらにどのような配列がエンリッチしているかを調べることでどのような転写因子がそこに結合するかを予測することができます. 例えばMed1のChIP-seqのピーク領域で結合モチーフ解析を行えばMed1が結合している場所にどのような転写因子が結合するかを予測することができます. 今回はこれをHOMERの`findMotifsGenome.pl`を使って行ってみます.  

これも下準備が必要でHOMERにゲノムをインストールする必要があります. 

```
perl $HOME/miniconda3/share/homer-4.10-0/configureHomer.pl -install mm10
```

miniconda3ディレクトリ以下の上のような位置にある`configureHomer.pl`でmm10をインストールします. 上記はホームディレクトリにminiconda3が入っている場合です. またhomer-4.10-0の部分は皆さんが`conda install homer`でHOMERをインストールした時期や環境で変わる可能性がありますので注意してください. （HOMERのバージョンが異なるかもしれないからです.）WSLの人でここでうまく行かなかった人は下のヘルプ2を参照してください.  

それでは`findMotifsGenome.pl`を走らせます.    

```
findMotifsGenome.pl Med1_data/Med1_TagDir/peaks.txt mm10 Med1_data/Med1_motif -size given -p 4
```

`findMotifsGenome.pl`の直後に`findPeaks`の結果のファイルを書きます. (他のピークコーラーの結果を使いたい場合は下のヘルプ3を参照してくだい.) その次にゲノムのバージョンを書きます. さらに結果をセーブするディレクトリの名前を書きます. `-size`オプションに入力した整数の長さの各ピークの中心を中心とする領域を対象に計算を行います. わかりにくいので例を挙げると`-size 100`とした場合は各ピークの中心から+/- 50 bpの範囲で計算を行います. 今回は`-size given`としてプログラムに渡したピーク領域それぞれの実際そのままの範囲で計算することにしました. `-p`には使用するスレッド数を入力します.  

計算が終わると`Med1_data`の中に`Med1_motif`というディレクトリができるので中を見てみましょう. 注目すべきは`knownResults.html`です. これをブラウザで開いてください. ここで計算が終わらずエラーが出た人も下のヘルプ2を参照してください.

![Motif.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/4db9bd02-a008-3315-f59e-4f09774edb1f.png)
 

このような画面になります. 今回の計算の場合, 自動的に用意されたランダムな領域に対して入力したピーク領域にHOMERのデータベース上のどのようなモチーフ配列がエンリッチしていたかをp値が小さい順に表示してくれています. ランダムな領域に対して計算を行うので結果は毎回ほんのすこしづつ変わります. この結果をみるとMed1のピーク領域にはKLF, OCT, SOXなどの幹細胞性維持に重要な転写因子のモチーフが多く存在し, Med1とこれらの転写因子が共局在していることを示唆しています.  

これとは別に`homerResults.html`では入力したピーク領域に多く存在していてHOMERのデータベース上のモチーフ配列と似ている配列についてまとめられています. 基本的には`knownResults.html`を確認すれば良いでしょう.  

`findMotifsGenome.pl`についてさらに詳しくは http://homer.ucsd.edu/homer/ngs/peakMotifs.html を御覧ください.

# ピーク領域同士の重なりを調べる

最後にピーク同士の重なりを調べる方法について紹介します. これにはHOMERの`mergePeaks`を使います. 今回はMed1 ChIP-seqのピーク領域とH3K27Ac ChIP-seqのピーク領域の重なりについて調べます.  

```
mergePeaks -d given Med1_data/Med1_TagDir/peaks.txt H3K27Ac_data/H3K27Ac_TagDir/regions.txt -prefix mergePeaks -venn venn.txt
```

`-d`を`given`にすると入力したピーク領域同士の重なりをそのまま素直に計算します. 入力するピーク領域を書き並べます. 3つ以上あってもよいです. `-prefix XXX`としておくとファイル名が`XXX`で始まる各ピークの重なっている領域を結合した領域と, 一つのピーク領域にしか存在しない領域を分けて出力してくれます. `-venn YYY.txt`としておくと`YYY.txt`という重なっているピーク領域の数, どちらか一方にしかないピーク領域の数をベン図を描くためにまとめた表を作ってくれます. オプションなどについて詳しくは http://homer.ucsd.edu/homer/ngs/mergePeaks.html を参照してください.  


このコマンドを実行すると`mergePeaks_H3K27Ac_data_H3K27Ac_TagDir_regions.txt`, `mergePeaks_Med1_data_Med1_TagDir_peaks.txt`, `mergePeaks_Med1_data_Med1_TagDir_peaks.txt_H3K27Ac_data_H3K27Ac_TagDir_regions.txt`, `venn.txt`というファイルが作られます. それぞれがH3K27Ac ChIP-seqのピークのみに存在する領域, Med1 ChIP-seqのピークのみに存在する領域, Med1 ChIP-seqとH3K27Ac ChIP-seqのピーク領域のうち重なっているものを結合した領域, ベン図を描くための表となっています.　　

Pythonのmatplotlibでベン図を描いてみましょう. matplotlib_vennというパッケージを使うので, 

```
conda install matplotlib-venn
```

してください. 次に何でもいいのでエディタを開いて以下のコードを書いてください.

```python
from matplotlib import pyplot as plt
from matplotlib_venn import venn2

#venn.txtを開いてMed1 ChIP-seqのみに存在するピーク数, 
#H3K27Ac ChIP-seqのみに存在するピーク数, 重なっているピーク数を確認してください.

venn2(subsets=(770, 25254, 2738), set_labels = ("Med1", "H3K27Ac"))
#subsets=(Med1のみ, H3K27Acのみ, 重なっているピーク)

plt.savefig("./venn.png")
```

書けたら`venn_plot.py`などの名前で保存して保存した場所で以下のコマンドを実行してください.  

```
python venn_plot.py
```

するとそのディレクトリに`venn.png`というファイルができるので開いてください.  

![venn.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/477619/f7617322-7215-9181-097d-feb213a0dd12.png)


転写を活性化する蛋白であるMed1のChIP-seqのピークの約8割が同じく転写活性マーカーであるH3K27AcのChIP-seqのピークと重なっています. にしてもH3K27Ac ChIP-seqのピーク数が多いですね. ぜひIGVで２つのデータを見比べてみてください.


# 終わりに

ここまで読んでくださってありがとうございます. ピークコールしてからの解析で今回ご紹介したものはChIP-seq解析のほんの一部です. みなさんはもうマッピング, ピークコールができるようになったのであとは目的に合わせた解析をHOMERやRやPythonなどを用いて行ってください. この記事を含めバイオインフォマティクスについての情報はいまやネットにあふれています. また今回と同様の手順でクロマチンの開閉状態についてゲノム全体で網羅的に調べるATAC-seq(Assay for Transposase-Accessible Chromatin Sequencing)の解析も可能です. この記事がみなさんの研究の助けになることを願っています. もし質問などございましたらできる限りお力になれればと思いますのでコメントお待ちしております！ 

# ヘルプ

## 1. `conda`コマンドが動かない！

おそらくminiconda3をインストールしたあとにパスが通っていないことが原因であると考えられます. 以下の操作を行ってください.  

```
cd #ホームディレクトリへ移動

vim .bash_profile #パスを記述する.bash_profileをvimで開く
```
開いたらまずescキーを押してください. つぎにiのキーを押してください. これでInsert modeで編集できるようになります. 以下を間違えないように入力してください.  

```
PATH=~/miniconda3/bin:$PATH
```

上記はminicondaがホームディレクトリにある場合です. 間違えていないかよく確認したらもう一度escキーを押してください. そして`:wq`と入力してenterを押してください.  

そのあとターミナルを再起動するか

```
source .bash_profile
```
を実行してください. これでパスが通って`conda`が動くはずです.

## 2. `configureHomer.pl -install mm10`がうまく行かない！  

おそらくWSLの人だけが起こりうるエラーなのですが`configureHomer.pl`を実行するのに必要なコマンドがインストールされていないことが原因かもしれません. 以下をすべて実行してください.

```
which gcc
which g++
which make
which perl
which zip
which gzip
which wget
```

このうち

```
/usr/bin/make
```

のようにそのコマンドへのパスが表示されなかったものがあれば

```
sudo apt install zip #zipが入っていなかったとき
```

を実行してインストールしてください. なかったものをすべてインストールしたあとでもう一度

```
perl $HOME/miniconda3/share/homer-4.10-0/configureHomer.pl -install mm10
```
を実行してみてください.

## 3. HOMERでHOMER以外のプログラムで作成したbedファイルを使う

HOMERでHOMER以外のプログラムで作成したbedファイルを使う場合は予めそれをHOMER bedファイルに変換しておくと安全です. それにはHOMERの`bed2pos.pl`を使用します.  

```
bed2pos.pl (変換したいbedファイル) > Converted_file.hb
```

HOMER bedファイルの拡張子は"hb"です.  

## 4. WSL でカレントディレクトリを Windows のエクスプローラーで開くには？  

```  
explorer.exe .  
```  

最後の `.` はカレントディレクトリという意味です.  

# References

- [Kagey MH *et al.,* Mediator and cohesin connect gene expression and chromatin architecture. Nature 467, 430–435 (2010).](https://www.nature.com/articles/nature09380)

- [Creyghton *et al.,* Histone H3K27ac separates active from poised enhancers and predicts developmental state. PNAS 107, 21931–21936 (2010).](https://www.pnas.org/content/107/50/21931.long)

- http://rnakato.hatenablog.jp/entry/2017/07/06/110926

- https://bi.biopapyrus.jp/rnaseq/qc/trimmomatic.html

- http://www.usadellab.org/cms/?page=trimmomatic

- https://bi.biopapyrus.jp/rnaseq/qc/fastqc.html

- https://bi.biopapyrus.jp/rnaseq/mapping/bowtie2/

# 免責事項  
本記事の正確性や妥当性につきまして，一切の保障はいたしません．また本記事を参考にして解析を行ったことで生じたあらゆる損害について，一切の責任を負いません．