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

今回は、全染色体を１行でマクロに可視化するGVコマンドを使います。なおGVはGlobal viewの略です。

parse2wig

今回はROADMAP web portalからダウンロードしたK562細胞のヒストン修飾データ一式を使います。以下のコマンドでtagAlignファイルをダウンロードします。（ダウンロードに時間がかかる場合は適宜サンプル数を減らして大丈夫です。）

$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-H3K4me3.tagAlign.gz
$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-H3K9me3.tagAlign.gz
$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-H3K27ac.tagAlign.gz
$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-H3K27me3.tagAlign.gz
$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-H3K36me3.tagAlign.gz
$ wget http://egg2.wustl.edu/roadmap/data/byFileType/alignments/consolidated/E123-Input.tagAlign.gz

ダウンロードしたファイルからparse2wigでWigデータを生成します。 ROADMAPのデータはgenome build hg19です。

for name in H3K27ac H3K27me3 H3K36me3 H3K4me3 H3K9me3 Input
do
    parse2wig -i E123-$name.tagAlign.gz -o E123-$name -gt genometable.txt -f TAGALIGN -binsize 100000 
done

# 上のfor文は以下のコマンドを実行したのと同じ
$ parse2wig -i E123-H3K27ac.tagAlign.gz -o E123-H3K27ac -gt genometable.txt -f TAGALIGN -binsize 100000
$ parse2wig -i E123-H3K27me3.tagAlign.gz -o E123-H3K27me3 -gt genometable.txt -f TAGALIGN -binsize 100000
$ parse2wig -i E123-H3K36me3.tagAlign.gz -o E123-H3K36me3 -gt genometable.txt -f TAGALIGN -binsize 100000
$ parse2wig -i E123-H3K4me3.tagAlign.gz -o E123-H3K4me3 -gt genometable.txt -f TAGALIGN -binsize 100000
$ parse2wig -i E123-H3K9me3.tagAlign.gz -o E123-H3K9me3 -gt genometable.txt -f TAGALIGN -binsize 100000
$ parse2wig -i E123-Input.tagAlign.gz -o E123-Input -gt genometable.txt -f TAGALIGN -binsize 100000

入力がtagAlign形式データのため、オプションに -f TAGALIGNを指定しています。また、GVコマンドはデフォルトbinsizeが100kbpのため、-binsize 100000 を指定し、100kbpのWigファイルを生成します。

DROMPAの実行

それでは可視化してみましょう。GVコマンドもPC_SHARPコマンドと指定方法はほぼ同様です。

$ drompa_draw GV \
$ -i parse2wigdir/E123-H3K4me3,parse2wigdir/E123-Input,H3K4me3 \
$ -i parse2wigdir/E123-H3K27ac,parse2wigdir/E123-Input,H3K27ac \
$ -i parse2wigdir/E123-H3K27me3,parse2wigdir/E123-Input,H3K27me3 \
$ -i parse2wigdir/E123-H3K36me3,parse2wigdir/E123-Input,H3K36me3 \
$ -i parse2wigdir/E123-H3K9me3,parse2wigdir/E123-Input,H3K9me3 \
$ -p K562 -gt genometable.txt

以下のような図が生成されます。 f:id:rnakato:20180216144039p:plain

これはChIP/Input比を100kbpの解像度で可視化したもので、ChIP/Input >1.0の領域はオレンジ、そうでない領域は灰色で表示されます。
ヘテロクロマチンマーカであるH3K9me3のS/N比は非常に弱いため、通常のピーク抽出ではほとんどピークが抽出されないこともしばしばありますが、このようなマクロな可視化で見ると、他のアクティブマーカであるヒストン修飾と排他的に局在するさまが視覚的によくわかります。染色体の免疫染色のようなイメージですね。

オプションをつければPC_SHARPと同様にリード数そのものを可視化することももちろん可能です。（以下の記事参照）

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

GC含量と遺伝子密度のグラフを追加

上の図だと少し殺風景なので、GC含量と遺伝子密度のグラフを追加してみましょう。 hg19のGC含量と遺伝子密度のデータはDROMPAのweb siteからダウンロード可能です。

$ wget http://www.iam.u-tokyo.ac.jp/chromosomeinformatics/rnakato/drompa/annotations/GC_hg19.tar.gz
$ wget http://www.iam.u-tokyo.ac.jp/chromosomeinformatics/rnakato/drompa/annotations/gene_density_hg19.zip
$ tar zxvf GC_hg19.tar.gz
$ unzip gene_density_hg19.zip

GC_hg19は500kbpごとのゲノム配列のGC含量、gene_density_hg19はその領域に含まれる遺伝子数です。これらをプロットするには以下のように指定します。

$ drompa_draw GV \
$ -i parse2wigdir/E123-H3K4me3,parse2wigdir/E123-Input,H3K4me3 \
$ -i parse2wigdir/E123-H3K27ac,parse2wigdir/E123-Input,H3K27ac \
$ -i parse2wigdir/E123-H3K27me3,parse2wigdir/E123-Input,H3K27me3 \
$ -i parse2wigdir/E123-H3K36me3,parse2wigdir/E123-Input,H3K36me3 \
$ -i parse2wigdir/E123-H3K9me3,parse2wigdir/E123-Input,H3K9me3 \
$ -p K562_withGCandGene -gt /home/Database/UCSC/hg19/genome_table \
$ -GC GC_hg19 -gcsize 500000 \
$ -GD genedensity -gdsize 500000

-GC, -GD オプションでそれぞれのデータ（ディレクトリ）を指定します。-gcsize, -gdsizeで、ウィンドウサイズを指定します。ここでは500kbpです。

結果は以下のようになります。

f:id:rnakato:20180216144049p:plain

H3K9me3が濃縮している領域は遺伝子数の少ないgene desertのような領域が多いことがざっと観察できます。 GC含量は、例えばリードの濃縮が極端にGC richに偏っているかどうか調べたい場合などに便利です。

ちなみに、EnsemblやUCSC genome browserで使われている↓のような簡単な染色体図も同時に可視化できるようにしたいのですが、著作権フリーな染色体模式図の画像ってどこかにあるんでしょうか・・・ご存知の方がおられましたら教えてください。

f:id:rnakato:20180216145832p:plain

補足：GC含量データ、遺伝子密度データの作成

GC含量、遺伝子密度のデータは自分で作成することもできます。ここではヒトゲノムhg19を例に解説しますが、他のゲノムの場合には適宜読み替えてください。

GC含量はDROMPA3のscriptフォルダに含まれる GCcount.plを使います。ここではGC_hg19フォルダを新たに作り、その中にデータを生成していくことにします。

$ mkdir GC_hg19  # フォルダ作成
$ winsize=500000  # ウィンドウサイズを500kbpに指定
$ for i in $(seq 1 22) X Y M # 1~22番染色体とX, Y染色体、ミトコンドリアを指定。
$ do
$    GCcount.pl chr$i.fa $winsize > GC_hg19/chr${i}_bs$winsize
$ done

GCcount.pl は単一のfastaファイルを入力とするので、染色体ごとに適用します。出力ファイル名は<染色体名>_bs<ウィンドウサイズ>である必要があり、染色体名はgenometable.txtと一致している必要があります。
このように作成し、-GC GC_hg19と指定すれば可視化可能です。任意のウィンドウサイズで可視化可能で、PC_SHARPコマンドでも可視化できます。

遺伝子密度データはDROMPA3のscriptフォルダに含まれる makegenedensity.pl を使います。実行にはgenome tableと、refFlat形式の遺伝子ファイルが必要です。持っていない場合はUCSCからダウンロードしましょう。

# 遺伝子ファイルのダウンロード　
$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/refFlat.txt.gz
$ gunzip refFlat.txt.gz  # 解凍
# 遺伝子密度データ作成
$ winsize=500000  # ウィンドウサイズを500kbpに指定
$ makegenedensity.pl genometable.txt refFlat.txt $winsize

入力にするrefFlatファイルをあらかじめ編集しておくことで、例えばコーディング遺伝子のみ・ノンコーディング遺伝子のみカウントするような使い方も可能です。