はじめに
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 SRAからシークエンスデータをダウンロードするのにつかうツール
- Trimmomatic アダプタートリミングをするツール
- fastQC fastqファイルのQCを行うツール
- Bowtie2 マッピングをするツール
- Picard PCR duplicateを除去するツール
- samtools 色々なところで使うツール
- deepTools データの補正に使うツール
- HOMER 色々できるツール
以下のように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) からシークエンスデータをダウンロードします. 今回はKagey MH et al.,(GSM560348)からマウスのES細胞のMed1のChIP-seqのデータとCreyghton et al.,(GSM594579)から同じくマウスのES細胞のH3K27AcのChIP-seqのデータを使用します. また, Med1と同じ論文からインプットのシークエンスデータ(GSM560357)も取得します. ()内の番号は各データのGEO Accession番号です. NCBIのGene Expression Omnibus(GEO)でこれらのデータのページ飛ぶのに必要な番号です.
さて, データをダウンロードする前にまずはこれらのデータのありかをブラウザ上で確認しておきましょう. Med1のデータを例に説明します. まずはGEOのサイトをブラウザで開いてください. 下の画像で赤く囲っている検索窓にGSM560348と打ち込んでください.
下の画像のようにこのデータについての色々な情報が書いてあるのでしっかり目をとおしてください.
では, このページの下の方にあるSRAと書かれたところの右にある番号をクリックしてください.
すると以下のようなページに行きます. 画像で赤く囲っている番号がSRR番号と言ってこのデータのアクセス番号でダウンロードする際に必要なものです.
論文中にはこれらの番号がどこかに書いてあってデータのありかがわかるようになっています.
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?
それではシークエンス結果のクリーニングを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のページの"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というファイルができるのでブラウザで開いてみましょう.
左のSummaryのところがほとんど緑のチェックになっていたらクオリティはOKです. このデータきれいすぎる... 各指標については https://bi.biopapyrus.jp/rnaseq/qc/fastqc.html で説明してくださっています. こちらのサイト, bioinformaticsさんには常日頃からお世話になっています. いつもありがとうございます.
残りの2つのデータについても同じオプションで行ってくださいね.
マッピング
いよいよBowtie2を使ってマッピングして行きます！と, その前にゲノムのインデックスをビルドしなくてはなりません. つまり, マッピングに必要なリファレンスゲノムの準備を行う必要があります.
UCSCのページから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言語の勉強をしておきましょう.（いいねよろしくお願いします！）
この記事では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のダウンロードページで自分のOSのインストーラーをダウンロードしてください. IGVはGUI（グラフィックが出てきてマウスとキーボードで操作するインターフェイス）ですのでWSLを使っているWindowsユーザーの方はここではWindows版を選んでください. ダウンロードしたインストーラーを起動して指示に従ってインストールしてください. するとデスクトップに以下のようなIGVのショートカットが作られます.
これをダブルクリックして起動させてください(起動に30秒ほどかかります). 起動後以下のようなウインドウが現れます.
今, hg19がロードされているのでmm10をダウンロードしてロードしましょう. 上の画面で赤く囲っているところの下矢印をクリックすると"More..."と出てくるのでそれをクリックしてください. そして"Mouse mm10"をクリックして左下の"Download Sequence"にチェックをいれて"OK"をクリック. これでmm10のダウンロードが始まります.
そろそろbamCoverageが終わったころですかね...？
終わったらMed1.bigwigをIGVのウィンドウにドラッグアンドドロップしてください.
上の画像のようにMed1 ChIP-seqのプロファイルが表示されたでしょうか? この状態ではすべての染色体を俯瞰して見ており, 詳細がわからないので赤く囲ってある検索窓に色々な遺伝子名前を入力してそのgene bodyの場所に飛んでみましょう. 以下に一つだけ例を示します.
Klf4は幹細胞性維持に重要な遺伝子ですがこのようにマウスのES細胞では転写を活性化する蛋白であるMed1が結合していることがわかります. こちらの画像では見やすいようにRefseq genesのトラックを引き上げています. また, 表示の縦軸のスケールがおかしい場合は"Med1.bigwig"という名前が書いてあるところで右クリックして"Autoscale"をクリックしてください.
残りの2つのデータについても同様にbigwigを作成してIGVで確認してみてください.
ピークコール
それではシグナルのピークを統計学的基準に基づいて検出するピークコールを行います. 今回はHOMERのfindPeaksを使います. これ以外にピークコーラーとしてよく使われているものにMACS2があります. 興味があれば結果を比べてみてください.
さて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で開くと見やすいです.
上の方に様々な情報が書かれている後で画像の下半分のようなピークの情報が書かれています. 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行目以降にさらに情報が入っている場合もあります. )
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にドラッグアンドドロップしてみてください.
このようにピークコールされた位置が表示されるはずです. わかりやすい...
結合モチーフ解析
転写因子はある程度決まった配列（結合モチーフ配列）に結合するのでピーク領域が与えられたとき, それらにどのような配列がエンリッチしているかを調べることでどのような転写因子がそこに結合するかを予測することができます. 例えば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を参照してください.
このような画面になります. 今回の計算の場合, 自動的に用意されたランダムな領域に対して入力したピーク領域に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
してください. 次に何でもいいのでエディタを開いて以下のコードを書いてください.
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というファイルができるので開いてください.
転写を活性化する蛋白である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
免責事項
本記事の正確性や妥当性につきまして，一切の保障はいたしません．また本記事を参考にして解析を行ったことで生じたあらゆる損害について，一切の責任を負いません．