Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

記事一覧

解析環境構築

環境変数PATHの通し方

pyenvでPython環境を構築する

バイオインフォマティクスのためのpython環境構築方法を考える

GitHubからプログラムをダウンロード・インストール

データ生成

 

常染色体と性染色体のみのゲノム配列ファイル genome.fa を作成する

2bit genome を作成する

genome tableを作成する

LiftOver: BEDファイルを異なるbuildへ変換

 

fastqデータ取得

SRAからfastqを取得する

ENA,DDBJからfastqを取得する

 

ゲノムマッピング

Readをゲノムにマッピング (その1) (2017/12/19 追記あり)

Readをゲノムにマッピング (その2)

Readをゲノムにマッピング (その3) 圧縮ファイルを入力にする方法

 

マップデータ操作

SAMtoolsとリダイレクト

SAMtoolsワンライナー覚書

 

ChIP-seq解析: DROMPA

DROMPA3: その1 インストール

DROMPA3: その2 parse2wig

DROMPA3: その3 ピーク抽出(peak calling)

DROMPA3: その4 マップリード分布の可視化その1

DROMPA3: その5 シェル変数を使う

DROMPA3: その6 ChIP/Input ratio 及び p値の可視化

DROMPA3: その7 -i オプション詳細

DROMPA3: その8 GVコマンドでのマクロな可視化

DROMPA3: その9 リードプロファイル

DROMPA3: その10 ヒートマップ

ChIP-seq解析: 品質評価

Library complexity (PCR bias)とは何か

S/N比の評価手法 その1

S/N比の評価手法 その2 Cross-correlation profile

S/N比の評価手法 その3 deepTools

S/N比の評価手法 その4 SSP

ChIP-seq解析: ピークを入力とする操作

BEDtoolsワンライナー覚書

2サンプル間ピーク比較

2サンプル間ピーク比較

2つのサンプルから得られたピークセットがどのくらい重なるのか調べたい!という時の方法です。

今回はBEDtoolsを使うやり方と、拙作のcompare_bsを使うやり方を紹介します。

ピークデータのダウンロード

ピークはBED形式であれば何でもよいのですが、ここでは例としてUCSCからH1細胞のCjunとMaxのピークファイルをダウンロードします。
(.narrowPeak という拡張子のファイルですが、これらはMACS2で得られたピークファイルで、BED形式として扱えます。)

 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak.gz
 $ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeAwgTfbsUniform/wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak.gz

解凍後、ピーク数を確認してみましょう。

 $  gunzip *.gz
 $  wc -l wgEncode*
   2148 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak
  11129 wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak
  13277 total

Cjunのピークはやや少ないようです。

bedtools

bedtoolsのintersectを使うと、共通するピーク領域を出力してくれます。その後数をカウントすればピーク数が判明します。

(Bedtoolsの詳細はこちら: BEDtoolsワンライナー覚書 - Palmsonntagmorgen

 $ bedtools intersect \ 
    -a wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak \
    -b wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak \
     > and.bed
 $  wc -l and.bed
377 and.bed

Cjunの2148ピークのうち、Maxと重なるものは377ピークになりました。

bedtools intersectは2つのピークが1塩基でも重なっている場合にoverlapと判定します。 また、サンプルAの1つのピークに大してサンプルBの3つのピークが重なる場合、出力されるピーク数は3つになります。コマンドに -wa オプションをつけると、1つのピークとして出力されます。
intersectで出力される領域については、以下に図解があります。
intersect — bedtools 2.27.0 documentation

compare_bs

compare_bsは以下の記事で説明したChIPseqTools の中にあります。

GitHubからプログラムをダウンロード・インストール - Palmsonntagmorgen

compare_bsを使うと以下のようになります。

 $ compare_bs -and \
    -1 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak \
    -2 wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak \
    > and.bed
 $ head and.bed
#file1: wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak
#file2: wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak
#num1: 2148     num2: 11129     num1_overlap: 373 (17.4%)       num1_notoverlap: 1775 (82.6%)   num2_overlap: 378 (3.4%)num2_notoverlap: 10751 (96.6%)
#peakwidth total1: 523816 bp    peakwidth total2: 3534874 bp    overlappeaks total: 70413 bp (13.44% / 1.99%)
chromosome      start   end     summit  length  enrich  intensity
chr10   80917404        80917527        80917465        123     0.00    149.81
chr22   39570815        39571059        39570937        244     0.00    142.53
chr18   9643685 9643754 9643719 69      0.00    133.02
chr6    20811494        20811738        20811616        244     0.00    132.72
chr1    205263766       205264010       205263888       244     0.00    116.91

compare_bsの出力では、冒頭にピーク数のスタッツが表示されます。 この場合、Cjunのピーク数2148、Maxのピーク数11129で、CjunのうちMaxと重なるものが373ピーク、MaxのうちCjunと重なるものが378ピークとわかります。
(CjunとMaxでそれぞれ重なったピーク数がずれているのは、ピークの重なりが1対1対応でないことに起因します。つまり1つのピークに対して複数の相手方ピークが重なる可能性があるということです。)

その下の行の #peakwidth total は、ピークの数でなく延べピーク幅で重なりをカウントしたものです。

compare_bs はデフォルトでは「サンプル1のピークのうち、サンプル2と重ならないもの」を出力し、-and オプションをつけると重なるものを出力します。 ですので、出力されるピークは bedtools intersect の -waオプションと同じく、サンプル1のピーク幅になります。

重なりの定義もbedtools と同じく、1塩基でも重なればoverlapと判定します。 -l 10000 とオプションをつけると、10kbpまで離れたピークをoverlapとして判定します。

 $  ~/git/ChIPseqTools/bin/compare_bs -and -l 10000  \
     -1 wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak  \
     -2 wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak \
     -nobs
#file1: wgEncodeAwgTfbsSydhH1hescCjunIggrabUniPk.narrowPeak
#file2: wgEncodeAwgTfbsSydhH1hescMaxUcdUniPk.narrowPeak
#num1: 2148     num2: 11129     num1_overlap: 689 (32.1%)       num1_notoverlap: 1459 (67.9%)   num2_overlap: 787 (7.1%)num2_notoverlap: 10342 (92.9%)
#peakwidth total1: 523816 bp    peakwidth total2: 3534874 bp    overlappeaks total: 6727767 bp (1284.38% / 190.33%)

ここでは-nobs オプションをつけて、ピークリストを出力せずに冒頭のスタッツのみを表示しています。 10kbpまで重なりの閾値をゆるめたことで、少し重なりが増えました。
(注:overlapped peakwidthは10kbp extensionも考慮してoverlapを出すので、パーセンテージは少しおかしくなってしまいます)

まとめ

bedtools, compare_bsで得られる結果を図にまとめると以下のようになります。

f:id:rnakato:20180806141533p:plain

得られるピークサイトはほぼ同じですが、compare_bsの有意点としては、statsを出してくれることとサンプル1と2のoverlap peak数をそれぞれ出してくれること、peak widthでもstatsを出してくれることがあります。
ピーク比較の時にはピークの幅は一定であると暗黙に仮定してしまいがちですが、幅の分布はかなり広がっているのが普通で、 1つの大きなピークが細切れになってカウントされていたりすることもあるので、ピーク数をカウントしてものを言いたい時には実際に得られたピークを確認する作業を忘れないようにしましょう。

GitHubからプログラムをダウンロード・インストール

NGS解析のための新規ツールは日々論文で発表されており、それらのほとんどは世界中の人が無償で利用可能なライセンス形態になっています。
今日はその中でも多くの人に利用されている「GitHub」に公開されたツールのインストール方法を紹介します。

オープンソースについて

今のBioinformatics界は極めてオープンソース志向が強く、「良いものや知見は共有されることでより良いものができる」という性善説的な発想に基づいています。bioRxivのようなopen preprintにも積極的です。
ので、新規のプログラムで論文を書く場合はそのプログラムがオープンソースで公開されていることが必須条件となっている論文誌がほとんどです。オープンにする際も、自分のラボページにzipファイルをアップするだけ、のようなあやしげな状態ではなく、GitHubBioconductorのような広く使われバージョン管理がしっかりしているコミュニティスペースに公開されることが推奨されます。Bioinformatics誌のレビュアーの一人は、「GitHubに公開されていないプログラムはその時点でrejectにする」とツイートしていました。

一方、たとえばNatureやCellのような大きな論文内で使わている個々の解析技術はin-houseで作られたスクリプトであることが多く、そういうのはあまり公開はされません。Wetのラボはソースの公開に消極的なところが多い印象です。

GitHub

このブログでは僕のツールを使うことが多いので、僕のアカウント(rnakato)を例に以下説明します。 僕のGitHubサイトはこちらです。DROMPAやSSPなど、今まで開発したツールやスクリプトを置いています。

github.com

GitHubについて

Git・GitHubとは非常に平たく言うと、作成したソースコードをネット上に公開して誰でもダウンロードできるようにするためのしくみです。ファイルのバージョン管理機能に優れており、複数の人が同時にソースを修正するような開発に特に向いています。
(どこかの新聞記事では「設計図共有サイト」とか言っててコーヒー吹いた
より詳しく知りたい場合は「GitHub 入門」でググってください。

GitHub上のツールをダウンロードするだけならばアカウント作成は大丈夫です。

GitとGitHubのちがい

Gitはファイルのバージョン管理のしくみそのもので、GitHubはGitをWeb上で利用するためのサービス(Gitホスティングサービス)です。 GitHub以外にもBitbucketやGitLabなどいくつかのGitホスティングサービスがあります。

Gitのインストール

gitとコマンドを打ってみて、「そんなツールないよ」とエラーが出たら、以下のコマンドでインストールしましょう。

 $ sudo apt install git

GitHubからインストール(スクリプトの場合)

スクリプトとは、そのファイルをダウンロードすればコンパイル無しでそのまま使えるソースコードです。 私が作っているスクリプトを以下で公開しています。(自分用なので、大半は使い物にはならないと思いますが、一部は有用と思います。)

GitHub - rnakato/script_rnakato

上のリンクを開くと↓のような画面になります。画面右にある"Clone or download"をクリックするとクローン用URLが表示されますので(図の赤枠部分)、コピーしてください(URL右のアイコンをクリックすると自動でコピーされます)。

f:id:rnakato:20180804142534j:plain

またはURL下の"Download zip"をクリックすると、圧縮ファイルとしてダウンロードすることもできます。

git clone <コピーしたURL> とコマンド実行すると、内容がクローン(ダウンロード)されます。

 $ git clone https://github.com/rnakato/script_rnakato.git
Cloning into 'script_rnakato'...
remote: Counting objects: 1674, done.
remote: Compressing objects: 100% (25/25), done.
remote: Total 1674 (delta 15), reused 17 (delta 7), pack-reused 1642
Receiving objects: 100% (1674/1674), 6.72 MiB | 3.53 MiB/s, done.
Resolving deltas: 100% (1002/1002), done.
Checking connectivity... done.
 $ ls  # ダウンロードされたフォルダを確認
script_rnakato

これで、script_rnakatoフォルダの中にあるスクリプトが全てダウンロードされました。

内容の更新

git cloneでクローンしたフォルダの内部でgit pull origin masterとコマンドを打つと、クローン後に加えられたGitHub上の更新内容を取り込むことができます。

 $ cd script_rnakato
 $ git pull origin master
From https://github.com/rnakato/script_rnakato
 * branch            master     -> FETCH_HEAD
Already up-to-date.

今はクローンしたばかりですので、当然修正はありません。

GitHubからインストール(プログラムの場合)

DROMPA3をダウンロード・インストールしてみます。 スクリプトと同じ要領ですが、C言語ソースのため、ダウンロード後にコンパイルが必要です。

GitHub - rnakato/DROMPA3: DROMPA version 3 (Peak-calling, Visualization, Normalization and QC for ChIP-seq analysis)

 $ git clone https://github.com/rnakato/DROMPA3
 $ cd DROMPA3
 $ make

GitHubからインストール(submoduleを含むプログラムの場合)

以下のChIPseqToolsはC++プログラムですが、レポジトリ内部にSSPのレポジトリを含んでいます。 これはsubmoduleという機能で入れ子構造になったGitHubレポジトリです。

GitHub - rnakato/ChIPseqTools

submoduleを含むレポジトリをクローンする場合は、--recursiveオプションをつけます。

 $ git clone --recursive https://github.com/rnakato/ChIPseqTools.git
 $ cd ChIPseqTools
 $ make

また、内容を更新したい場合は以下のようにします。

 $ git pull origin master    # ChIPseqToolsレポジトリの更新
 $ git submodule foreach git pull origin master # ChIPseqTools内のSSPレポジトリの更新

まとめ

3種類のパターンに分けてGitHubからのツールのインストール方法を紹介しました。インストール方法は通常はページ下部のREADMEに書いてあるので、それを見れば大丈夫です。 NGS解析のツールの大半はGitHub上で公開されているので、GitHubを使いこなせるようになると自由度が各段にアップします。

またGitHub上に自分のアカウントを作成すれば、自分のスクリプトをアップすることもできるようになります。 GitHubは常に最新バージョンが保てる点がメリットなので、例えば職場と自宅でスクリプトを共有したいというような目的にも使えますし、 外部の人にスクリプトを送るのにも便利です。
(Dropboxでも共有はできますが、ファイル喪失の危険はありますし)

なおGitHubは無料プランだとファイル公開しか選べません。もしWeb上にアップしたいけれども公開にはしたくないという場合にはBitbucketかGitLabを利用しましょう。

LiftOver: BEDファイルを異なるbuildへ変換

公開されているゲノム配列は現在も更新中であるため、いくつかのバージョン (build) があります。 humanだとhg18, hg19, hg38などがあり、hg38が現時点で最新です。
NGS解析をするうえでは全ての解析データのbuildを統一する必要がありますが、「既存論文のデータの一部を用いたい」というような場合、buildが合わずに困ることがあります。 「自分はhg38で統一しているのに、既存論文のBEDファイルはhg19しか提供されてない」というような場合ですね。

そういった場合には LiftOver というツールを使うと、バージョンに合わせてBEDファイルを再生成してくれます。 詳細は以下のリンクを参照してください。

LiftOver - Genome Analysis Wiki

(話がそれますが、hg38がリリースされて既に久しいのですが、2018年の論文でも普通にhg19が使われています。 おそらくhg19で生成されたデータベース、アノテーションの類があまりに多すぎるため、 hg38にアップデートできない状態なのではないかと思います。ROADMAPもJuicerToolsもhg19ですしね。 既存データベースのデータを多く使う実験を組む際にはbuildをhg19にするのも手かもしれません。)

LiftOver (web)

hg19のBEDファイル(sample.bed)をhg38に変換する例です。
最も簡単なのは、UCSC web browser上で変換する方法です。

Lift Genome Annotations

デフォルトではhg19 -> hg38への変換になっていますので、そのまま使いましょう。 sample.bedをアップロードするか、bedの内容を直接テキストボックスに貼りつけてsubmitします。

submit すると、下にResultsができ、生成されたBEDファイルをダウンロードすることができます。

LiftOver (コマンドライン)

LiftOverはコマンドラインツールでも提供されています。 変換したいファイルが大量にあるような場合はこちらのほうが便利でしょう。

以下のようにしてバイナリファイルをダウンロードしましょう。

wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver
chmod +x liftOver  # 実行権限の付与

変換にはLiftOverツールそのものの他に、変換用のchainfileが必要になります。

以下のディレクトリから、hg19 -> hg38 用のchainfile hg19ToHg38.over.chain.gz をダウンロードしましょう。

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/liftOver/

ダウンロード後、解凍しておきます。

gunzip hg19ToHg38.over.chain.gz 

liftOver コマンドは、 liftOver <入力ファイル> <chainfile> <出力ファイル> <変換失敗したサイトを格納するファイル>というかたちで実行します。 今回の場合は以下のコマンドで変換完了です。

liftOver sample.bed hg19ToHg38.over.chain output.bed unlifted.bed

unlifted.bedには、何らかの理由で変換できなかったbedサイトが格納されます。変換できない理由は以下を参照してください。

https://genome.sph.umich.edu/wiki/LiftOver#Various_reasons_that_lift_over_could_fail

注意点

LiftOverを使っての変換は同一生物種の異なるbuildが想定されており、異なる生物種への変換は可能ではあるものの推奨はされていません。

また、変換失敗するサイトが発生することからもわかるように、変換は厳密なものではありません。 例えば、hg19でピークコールしたものをhg38に変換したピークリストと、hg38でピークコールしたピークリストはけっこう違います。遺伝子のアノテーションなどもかなり変わっています。
ので、大規模に統合したい場合には面倒でもfastqファイルから(条件を統一して)解析しなおす方がおすすめです。open chromatin regionのリストが欲しい、くらいであれば問題ないとは思います。

DROMPA3: その10 ヒートマップ

今回はDROMPAを用いたヒートマップの描画について説明します。

HEATMAPコマンド(TSS周辺)

前回の記事では、DROMPAを用いたリードプロファイルの描画について解説しました。

rnakato.hatenablog.jp

使ったコマンドは以下のようなものです。(詳しくは前回の記事を参照してください)

drompa_draw PROFILE \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562 -gene refFlat.txt -gt genometable.txt 

ヒートマップの描画はプロファイルと非常に似ています。 上のコマンドを以下のように少し変えるだけでヒートマップのコマンドになります。

drompa_draw HEATMAP \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p K562 -gene refFlat.txt -gt genometable.txt -png

変えた部分は、PROFILEをHEATMAPに変更したことと、ファイル名(-pオプション)、-pngオプションを付け加えたことです。
デフォルトではpdfでファイルが出力されますが、ヒートマップをpdfで出力するとしばしばファイルサイズが大きくなりすぎる&そこまで高解像度は必要ないため、-pngオプションを付けて.pngで出力させた方が扱いやすくなります。

PROFILEコマンドと同様、デフォルトでは入力された遺伝子の転写開始点周辺5kbを可視化します。 出力されたK562.heatmap.pngを開いてみましょう。 f:id:rnakato:20180709155123j:plain:w350

全遺伝子をヒートマップ化しているのでかなり縦長になります。ここに表示したのは一部です。
H3K4me3はTSSまわりに濃縮が見られる一方、H3K27me3, H3K36me3は強い濃縮が見られないことがわかります。

しかし、色のスケールを変えることで違いが出るかもしれません。 -scale_tag 5 オプションを追加して、色のMAX値を5に変更してみましょう。

drompa_draw HEATMAP \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p K562-2 -gene refFlat.txt -gt genometable.txt -png -scale_tag 5

f:id:rnakato:20180709160156j:plain:w350

どうやら、H3K27me3よりはH3K36me3 の方がTSSまわりに濃縮があるようです。

HEATMAPコマンド(ピーク周辺、ソートなし)

HEATMAPコマンドもPROFILEコマンド同様、ピーク周辺を可視化することができます。 また、複数のピークファイルを入力とすることができます。

ここではENCODEのオープンクロマチン領域を利用してみます。まずデータをダウンロードしましょう。

wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeOpenChromSynth/wgEncodeOpenChromSynthK562Pk.bed.gz
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeOpenChromSynth/wgEncodeOpenChromSynthHepg2Pk.bed.gz
# 解凍
gunzip wgEncodeOpenChromSynthK562Pk.bed.gz
gunzip wgEncodeOpenChromSynthHepg2Pk.bed.gz
# ピーク数が多いので、最初の1000行を抽出
head -n 1000 wgEncodeOpenChromSynthHepg2Pk.bed > wgEncodeOpenChromSynthHepg2Pk.1000.bed
head -n 1000 wgEncodeOpenChromSynthHepg2Pk.bed > wgEncodeOpenChromSynthHepg2Pk.1000.bed

DROMPAコマンドは以下のようになります。 -bed コマンドは、ファイル名とラベル(出力図y軸に表示されるもの)をコンマで区切って指定します。

drompa_draw HEATMAP \
            -i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
            -i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
            -i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
            -bed wgEncodeOpenChromSynthK562Pk.1000.bed,open_K562 \
            -bed wgEncodeOpenChromSynthHepg2Pk.1000.bed,open_HepG2 \
            -p K562-3 -ptype 4 -gt genometable.txt -png

f:id:rnakato:20180709163118p:plain:w340

2つのbedファイルで指定された領域が線で区切られて出力されました。

HEATMAPコマンド(ピーク周辺、ソートあり)

y軸はデフォルトでは入力ファイルの通りの順で出力されますが、 以下のように-hmsort オプションをつけると、ピーク強度でソートされます。

drompa_draw HEATMAP \
            -i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
            -i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
            -i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
            -bed wgEncodeOpenChromSynthK562Pk.1000.bed,open_K562 \
            -bed wgEncodeOpenChromSynthHepg2Pk.1000.bed,open_HepG2 \
            -p K562-4 -ptype 4 -gt genometable.txt -png -hmsort 1

f:id:rnakato:20180709163005p:plain

-hmsort 1とすると、1番目のサンプル(ここではH3K4me3)のピーク強度に基づいてソートされます。また、中心のリード数が0のサイトは除外されます。
こうしてみると、K562のオープンクロマチンサイトだけでなく、HepG2のオープンクロマチンサイトにもH3K4me3の濃縮が見られます。両方の細胞で共通するオープンクロマチンサイトなのかもしれません。

ここではH3K4me3にしか強い濃縮が見られませんが、共起するピークとそうでないピークに分けてベン図と組み合わせてヒートマップを使うと効果的にピークの濃縮を可視化することができます。

ヒートマップの特徴

PROFILEコマンドでは平均値(と信頼区間)を可視化しますが、この場合、全入力領域のうち全てに濃縮が出ているのか、あるいは半々程度にばらついているのか区別ができません。ヒートマップを使うと、全体のどのくらいの割合に濃縮が出ているのか、サンプル間で共通性・排他性があるのかなどを視覚的にとらえることができます。
一方、ヒートマップは統計的な処理をしているわけではないので、通常のリードプロファイルと同様、特徴的な領域を直感的に例示する場合などには有効ですが、全ゲノムを定量的に分析するような用途には使えません。

余談ですが、ヒートマップの描画はpythonの方がきれいにできるので、pythonに修正しようか考え中です。。

DROMPA3: その9 リードプロファイル

DROMPA3を用いたリードプロファイルの描画です。
ここでプロファイルと呼んでいるものは、遺伝子回り、あるいはピーク回りにおけるマップリードの平均分布のことで、aggregation plotと呼ばれることもあります。
全遺伝子や全ピークの平均値として見ることにより、通常のピーク抽出では得られない微小なリードの濃縮を捉えたり、発現遺伝子群と非発現遺伝子群の間でのリード濃縮の違いを可視化することができます。

インストール

DROMPA3のインストール方法はこの記事を参照してください。

PROFILEコマンドは内部でRスクリプトを実行しますので、Rが必要になります。

Genome table作成

DROMPA3の実行にはGenome tableファイルが必要になります。 Genome tableの作成は以下の記事を参照してください。

genome tableを作成する - Palmsonntagmorgen

遺伝子アノテーション

可視化する際に必要になる遺伝子アノテーションデータ(refFlat形式)をUCSCのサイトからダウンロードします。

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/refFlat.txt.gz
$ gunzip refFlat.txt.gz  # 解凍

parse2wig

今回もENCODE projectにあるK562のヒストン修飾群、遺伝子アノテーションを用います。

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwHistone/wgEncodeUwHistoneK562H3k4me3StdAlnRep1.bam
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwHistone/wgEncodeUwHistoneK562H3k27me3StdAlnRep1.bam
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwHistone/wgEncodeUwHistoneK562H3k36me3StdAlnRep1.bam
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwHistone/wgEncodeUwHistoneK562InputStdAlnRep1.bam
$ parse2wig -i wgEncodeUwHistoneK562H3k4me3StdAlnRep1.bam -o H3K4me3 -gt genometable.txt -f BAM
$ parse2wig -i wgEncodeUwHistoneK562H3k27me3StdAlnRep1.bam -o H3K27me3 -gt genometable.txt -f BAM
$ parse2wig -i wgEncodeUwHistoneK562H3k36me3StdAlnRep1.bam -o H3K36me3 -gt genometable.txt -f BAM
$ parse2wig -i wgEncodeUwHistoneK562InputStdAlnRep1.bam -o Input -gt genometable.txt -f BAM

PROFILEコマンド

drompa_draw PROFILE とタイプするとヘルプが表示されます。

$ drompa_draw PROFILE
Usage: drompa_draw PROFILE [options] -p <output> -gt <genometable> -i <ChIP>,<Input>,<name> [-i <ChIP>,<Input>,<name> ...]
       <output>: Name of output files
       <genometable>: Tab-delimited file describing the name and length of each chromosome
       -i: specify ChIP data, Input data and name of ChIP sample (separated by ',', <Input> and <name> can be omitted)

Options:
       (以下オプションが続く)

まずはデフォルトで実行してみましょう。コマンドは以下のようになります。

drompa_draw PROFILE \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562 -gene refFlat.txt -gt genometable.txt 

このコマンドにより以下のprofile.K562.pdfが生成されます。

f:id:rnakato:20180524223541p:plain:w350

このプロファイルは遺伝子の転写開始点(TSS)周辺5kbp (±2.5kbp) の平均マップリード数を可視化したものです。ラインが各サンプルの平均リード分布、陰になっている領域が95%信頼区間を表しています。
この図を見ると、H3K4me3はTSS周辺にリードの濃縮がみられる一方、H3K27me3とH3K36me3は濃縮がみられないことがわかります。

オプション

-cwオプションをつけると可視化する周辺長の長さを変えることができます。±5kbpのプロファイルを生成するには以下のようにします。

drompa_draw PROFILE \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562 -gene refFlat.txt -gt genometable.txt -cw 5000

f:id:rnakato:20180524223642p:plain:w350

デフォルトではChIPサンプルのリードの平均値を描画しますが、-stype 1 をつけるとChIP/InputのEnrichmentの平均値を描画します。

drompa_draw PROFILE \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562 -gene refFlat.txt -gt genometable.txt -stype 1

f:id:rnakato:20180524223834p:plain:w350

この例だとChIP readの時との差がわかりにくいですが、broad markや酵母の複製解析などでは差が顕著になることもあります。
(これは各TSSのChIP/Input の平均値であり、ChIP平均/Input平均のEnrichmentではないことに注意。)

遺伝子全長プロファイル

次に、遺伝子全長を100分割したプロファイルを作ってみましょう。 -ptype 3をつけて実行します。

drompa_draw PROFILE -ptype 3 \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562.ptype3 -gene refFlat.txt -gt genometable.txt 

f:id:rnakato:20180519134337p:plain:w400

このプロファイルは各遺伝子について長さを一定に正規化したうえで遺伝子長を100分割し、遺伝子内(gene body)でどのように濃縮しているかを見たものです。 x軸の0がTSS、100がTES(転写終了点)で、その両側に遺伝子長と同じ長さだけ伸長した領域について計測します。 (従って遺伝子ごとに計測している領域の幅が異なることに注意)
100分割した1つのウィンドウがwigファイルの1bin (100bp)より大きい場合、その領域内の全binの平均値を出力します。

TSS周辺の可視化ではほとんど平らだったH3K36me3ですが、遺伝子全体で見るとgene bodyに広くゆるやかな濃縮が広がっていることがわかります。H3K36me3は転写中のgene bodyに入る修飾ですので、これは妥当な結果と言えます。
H3K27me3の方はTSS周辺に微小な濃縮があるようですが、こちらは意味があるかどうか微妙なところです。

領域内でのリード数正規化

3つのサンプル間で縦軸がずれている(平均リード値が異なる)のが気になる場合は、-ntype 1オプションを追加すると、観測している領域内の総リード数でサンプルごとに正規化します。

drompa_draw PROFILE -ptype 3 -ntype 1 \
-i parse2wigdir/H3K4me3,parse2wigdir/Input,H3K4me3 \
-i parse2wigdir/H3K27me3,parse2wigdir/Input,H3K27me3 \
-i parse2wigdir/H3K36me3,parse2wigdir/Input,H3K36me3 \
-p profile.K562.ptype3 -gene refFlat.txt -gt genometable.txt 

f:id:rnakato:20180519151326p:plain:w350

完全ではありませんが、遺伝子外領域の高さがだいたい同じになりました。

まとめ

遺伝子全長の可視化は、H3K36me3のようにgene bodyに入るヒストン修飾や、RNA pol2のようにgene bodyに広く分布する(伸長する)タンパク質のChIP-seqを可視化するのに有効です。 今回はTSS周辺と遺伝子全長について可視化しましたが、-ptype 2でTES周辺、-ptype 4で指定したピークリスト周辺の濃縮を可視化することができます。

また、同時に生成されるRスクリプトprofile.K562.R を編集の後以下のように実行することで、ラベルやy軸のスケールを変更したpdfファイルを再生成することも可能です。

 R --vanilla < profile.K562.R

プロファイルの場合全体の何割で濃縮が入っているかはわかりませんが、信頼区間の幅を見ることである程度推測することが可能です。

S/N比の評価手法 その4 SSP

時間がかかってしまいましたが、やっとSSPの登場です。

この記事は以下の記事の続きです。

S/N比の評価手法 その1 - Palmsonntagmorgen

S/N比の評価手法 その2 Cross-correlation profile - Palmsonntagmorgen

S/N比の評価手法 その3 deepTools - Palmsonntagmorgen

SSPとは

SSP (Strand-shift profile)は Cross-correlation profile (上記その2参照)を拡張したもので、相関係数を取るかわりにJaccard indexを使って順鎖と逆鎖のマップリード分布の一致度を計算します。

Jaccard indexを用いるメリットとして以下が挙げられます。

  • 相関係数に比べて計算が高速。
  • PCR biasをフィルタした状態ではマップリード数の値が0か1しかないため、相関係数は不適。
  • unmap-unmapとmap-mapを両方評価する相関係数に対し、Jaccard indexはmap-mapペアのみを評価する。

SSPもいくつか機能があるのですが、今回はS/N比の計測機能のみに注目します。

SSPインストール

aptで必要なライブラリ群をインストールした後、ソースコードGitHubからダウンロードの後コンパイルするとsspの実行ファイルが作成されます。

$ sudo apt install git build-essential libboost-all-dev \
   libgsl-dev libz-dev samtools
$ git clone https://github.com/rnakato/SSP.git
$ cd SSP
$ make

以下のコマンドで実行ファイルにパスを通します。
(SSPをダウンロードしたディレクトリが $HOME/my_chipseq_exp/SSP の場合)
参考: 環境変数PATHの通し方 - Palmsonntagmorgen

$ export PATH = $PATH:$HOME/my_chipseq_exp/SSP/bin

sspとタイプしてヘルプが表示されればインストール成功です。

$ ssp                                                                                                             11:30:18

SSP v1.1.1
===============

Usage: ssp [option] -i <inputfile> -o <output> --gt <genome_table>
Use --help option for more information on the other options

Strand-shift profileの描画

その2 CCP の時と同じCTCFとInputのサンプルを使って、Strand-shift profileを描画してみましょう。

$ bam=wgEncodeUwTfbsHelas3CtcfStdAlnRep1.bam
$ ssp -i $bam -o CTCF -p 4 --gt genometable.hg19.txt --mptable mptable.hg19.txt
$ bam=wgEncodeUwTfbsHelas3InputStdAlnRep1.bam
$ ssp -i $bam -o Input -p 4 --gt genometable.hg19.txt --mptable mptable.hg19.txt

-i で入力サンプルを指定します。BAM, SAM, tagAlign, Bowtieフォーマットが入力可能で、ファイルのフォーマットはファイル名から自動で判定されます。
--gt はgenome tableファイル、 --mp はマップ可能な(mappable) ゲノム領域の長さを記したファイルです。いずれもインストールされたSSPフォルダのdata ディレクトリの中にあります。-pは使用するCPU数です。

出力ファイルはsspoutディレクトリの下に生成されます。CTCF.jaccard.pdf, Input.jaccard.pdf を見ると以下のようになっています。

CTCF.jaccard.pdf f:id:rnakato:20180427132049p:plain

Input.jaccard.pdf f:id:rnakato:20180427132126p:plain

左図はCCPとほぼ同じスケールですが、CCPと同じようなプロファイルが得られていることがわかります。 右図はx軸(順鎖と逆鎖のずれ)のスケールを1Mbpまで伸ばした時のプロファイル(log scale)です。
(CCPでは計算量の問題でx軸が1kbp程度、5bpおきの計測となっていますが、SSPは計算が高速なためsingle bp resolutionで遠くまで計算することができます)

左下にスコアが出ていますが、NSCS/N比の値として用いられるスコアです。CTCFは42.9, Inputは1.42 となっており, S/N比の差を正しく反映していそうです。

SSP, CCP, deepToolsの比較評価

それでは、SSP, CCP, deepToolsから得られるS/N比のスコアを比較してみましょう。
ここでの実験では以下のようにします。

  • sharp peak, broad peak, input の三種類のサンプルを用いる。
  • 1000万リード、2000万リードをランダムに抽出したファイルをそれぞれ作り、スコアが同じになるか比較(read depthの影響)。

サンプル数が多くなりますので、シェル変数やfor文を利用していきます。シェルスクリプトに不慣れな方は以下の記事も参考にしてください。

DROMPA3: その5 シェル変数を使う - Palmsonntagmorgen

データダウンロード

データはROADMAP projectのE018(iPS-15b cells)から、H3K4me3 (sharp), H3K9me3 (broad), Inputをダウンロードします。

for prefix in E018-H3K4me3 E018-H3K9me3 E018-Input; do
    wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/$prefix.tagAlign.gz
done

データ変換

ダウンロードされたサンプルは3000万マップリードありますが、ここから1000万、2000万リードを抽出したファイルを生成します。生成には SSP/scripts フォルダにある randompickup4bed.pl を使います。

for prefix in E018-H3K4me3 E018-H3K9me3 E018-Input; do
    file=$prefix.tagAlign.gz
    gunzip $file    #  gzファイルを解凍 $prefix.tagAlign ができる
    for num in 10 20; do
         # $prefix.tagAlign から ${num}000000 リードをランダムに抽出
         randompickup4bed.pl $prefix.tagAlign ${num}000000 > $prefix.${num}mil.tagAlign  
         # $prefix.${num}mil.tagAlign を圧縮 $prefix.${num}mil.tagAlign.gzができる
         gzip $prefix.${num}mil.tagAlign
    done
done

deepToolsは sorted bam 形式しか受け付けませんので、tagAlignファイルをsorted bamに変換します。

for prefix in E018-H3K4me3 E018-H3K9me3 E018-Input; do
    for num in 10 20; do
       bedToBam -i $prefix.${num}mil.tagAlign.gz -g genometable.txt | samtools sort -@4 > $prefix.${num}mil.sort.bam
       samtools index $prefix.${num}mil.sort.bam
    done
done

完了すると、以下のようにファイルができているはずです。

$ ls
E018-H3K4me3.10mil.sort.bam      E018-H3K9me3.10mil.sort.bam      E018-Input.10mil.sort.bam
E018-H3K4me3.10mil.sort.bam.bai  E018-H3K9me3.10mil.sort.bam.bai  E018-Input.10mil.sort.bam.bai
E018-H3K4me3.10mil.tagAlign.gz   E018-H3K9me3.10mil.tagAlign.gz   E018-Input.10mil.tagAlign.gz
E018-H3K4me3.20mil.sort.bam      E018-H3K9me3.20mil.sort.bam      E018-Input.20mil.sort.bam
E018-H3K4me3.20mil.sort.bam.bai  E018-H3K9me3.20mil.sort.bam.bai  E018-Input.20mil.sort.bam.bai
E018-H3K4me3.20mil.tagAlign.gz   E018-H3K9me3.20mil.tagAlign.gz   E018-Input.20mil.tagAlign.gz
E018-H3K4me3.tagAlign            E018-H3K9me3.tagAlign            E018-Input.tagAlign

実行

生成されたファイルを使ってSSP, CCP, deepToolsを実行します。 かなり時間がかかりますので、気長に待ちましょう。

for prefix in E018-H3K4me3 E018-H3K9me3 E018-Input; do
    for num in 10 20; do
           head=$prefix.${num}mil
           # SSP
           ssp -i $head.tagAlign.gz -o SSP-$head --gt genometable.txt -p 4 --mptable mptable.txt

           # CCP
           Rscript run_spp.R -p=4 -c=$head.tagAlign.gz -i=$head.tagAlign.gz \
           -savn=CCP-$head.narrowPeak -savr=CCP-$head.regionPeak \
           -out=CCP-$head.resultfile -savp=CCP-$head.pdf

           # deepTools
           plotFingerprint -b $head.sort.bam --JSDsample $head.sort.bam --ignoreDuplicates -p 4 \
           -plot DT-fingerprint.$head.png --outQualityMetrics DT-$head.txt
    done
done

結果表示

結果を見てみましょう。ひとつひとつファイルを確認するのは時間がかかるので、シェルスクリプトでぱぱっとやってしまいます。 (それぞれの行が何をしているのか考えてみてください)

echo -e "Sample\tSSP\tCCP\tdeepTools"
for prefix in E018-H3K4me3 E018-H3K9me3 E018-Input; do
    for num in 10 20; do
        head=$prefix.${num}mil
        echo -en "`tail -n1 sspout/SSP-$head.stats.txt | cut -f1,6`\t"
        echo -en "`cut -f9 CCP-$head.resultfile`\t"
        echo -e "`tail -n1 DT-$head.txt| cut -f9`"
    done
done

表示される結果は以下のようになります。

Sample SSP CCP deepTools
E018-H3K4me3.10mil  13.1502 1.429187    0.42114196480846744 
E018-H3K4me3.20mil  13.332  1.427474    0.4537138959592519  
E018-H3K9me3.10mil  3.37982 1.062391    0.2666993695154958  
E018-H3K9me3.20mil  3.36639 1.062076    0.29083049547676265 
E018-Input.10mil    1.31723 1.012775    0.1975224100365595  
E018-Input.20mil    1.29461 1.012816    0.20723845047525988 

タブ区切りの表になっているのですが、列の表示がずれていたり、有効桁数表示が多かったりして見づらいので(やたら桁数表示が多いツールってセンスがないですよね!)、エクセルに貼りつけて棒グラフにしてみましょう。

ぱっ!

f:id:rnakato:20180427145753p:plain:w400

f:id:rnakato:20180427145807p:plain:w400

f:id:rnakato:20180427145826p:plain:w400

ポイントをまとめると、

  • SSPはH3K9me3, H3K4me3ともInputより高いスコアを出しています。リード数にも依存していません。
  • これに対してCCPはH3K9me3とInputの差がほとんどありません。これが「sharp peakにしか使えない」と言われる所以です。
  • deepToolsはSSPと同じくH3K9me3, H3K4me3がInputより高いスコアを出していますが、リード数に依存してスコアが変わってしまっており、スコアがリード数に依存していることがわかります。
    • (2つしかないとわかりにくいですが、1000万リードよりも少なくなると急激に値が減少します)。

まとめ

以上、簡単ですが、SSPが他のツールよりも優れている点を比較評価してみました。 今回は3サンプルのみの簡単な実験ですが、論文では1000以上のChIP-seqサンプルを使った評価をしているので、興味があれば参照してください。

Sensitive and robust assessment of ChIP-seq read distribution using a strand-shift profile | Bioinformatics | Oxford Academic