bigwigの作成と表示

JBrowseでbamを表示させるとリードがずらーっと並ぶ
実際にいくつあるのかわかりやすくするために
ヒストグラムで表示したい

bamをbigwig形式にする
bam2wigとwigToBigWigを使う
ここを参考にした
bamからbigWigとWiggle Formatに変換するツール - macでインフォマティクス

コマンドはここからダウンロード
bam2wig: bamからwigを作成
GitHub - MikeAxtell/bam2wig: Conversion of a BAM alignment to wiggle and bigwig coverage files, with flexible reporting options
wigToBigWig: wigからbigwigを作成
Index of /admin/exe/linux.x86_64
両方にパスを通す
wigToBigWigにパスが通っていたら
bam2wigは自動でwigToBigWigも実行する

変換はbam2wigにbamを読み込ませるだけ

$ bam2wig hogehoge.bam

こんな名前のディレクトリhogehoge_bam2wigが
作られてbigwigも出力される

できたbigwigファイルをjbrowseで表示するには
jbrowseのデータのあるディレクトリにコピーして
スクリプトで取り込ませる

How do I set up a BigWig file?
JBrowse FAQ - GMOD

$ add-bw-track.pl --label mybw --bam_url file.bw

Example BigWig-based Wiggle XY-Plot Track Configuration
JBrowse Configuration Guide - GMOD

例えば作ったファイルhogehoge.bigwigが
以下のディレクトリに存在する場合

$ /var/www/html/jbrowse/hogehoge/bw

コマンドを以下のように実行した
trackList.jsonのあるディレクトリに移動して
スクリプトを実行するのが良さそう

$ cd /var/www/html/jbrowse/hogehoge/
$ ../bin/add-bw-track.pl --label hogehoge_bigwig --bw_url bw/hogehoge.bigwig --in ./trackList.json

trackList.jsonをvimで編集
hogehoge.bigwigのあるエントリーで
typeの項目を書き換える
"type" : "JBrowse/View/Track/Wiggle/XYPlot"

登録した結果はこんな感じ

青のヒストグラムがbigwig
ヒストグラムの下の部分は
ヒトゲノムのgff3を表示している
ちなみに異なる3種類のbigwigファイルを表示しているので
ヒストグラムが3段ある

ヒトゲノムのfastaファイルとgffは
以前作ったファイルとは変わっている
前にもちょっと書いたけど公開されているヒトゲノムは
ヘッダーの書き方に違いがある
GRCh38とhg38の違い(含むミトコンドリア)
これはhisatで作ったbamファイルのラベルが変わってくる

あと、gtfファイルはUCSCのゲノムブラウザを使って作成する
Gene annotation データを用意する（gtf形式） - Palmsonntagmorgen
Table Browser
gtfファイルはgffreadを使ってgffに変換する

$ gffread -E hg38_ucsc.gtf -o- > hg38_ucsc.gff3

ちなみに逆もできる

$ gffread rename_UCSC.hg38.gff3 -T -o- > rename_UCSC.hg38.gtf

このgffファイルにはシークエンスデータのパッチも含まれる
これが意外に厄介でパッチで追加された配列に重複があると
いろいろ変換ができなくなる

問題ありそうだけど重複部分を削除することにした
良くないんだろうけど。。。
ゲノムデータのパッチの取扱がよくわからない。

重複部分の削除とJBrowseの登録はこんな感じ

$ cat UCSC.hg38.gff3 | awk '(a[$9]++ < 1000) && ($1 !~/[_alt]/) {print}' > rename_UCSC.hg38.gff3
$ bin/flatfile-to-json.pl --gff ~/usr/data/db/rename_UCSC.hg38.gff3 --trackType CanvasFeatures --out hogehoge --trackLabel ref_hg38_ch

詳しくはまたあとで。

20 April 2019追記
作り直したhg38のgtfファイルについてはここにまとめてある
mecobalamin.hatenablog.com
このファイルをJBrowseに登録したり
マッピングに使ったりした

追記ここまで