ラベル 視覚化 の投稿を表示しています。 すべての投稿を表示
ラベル 視覚化 の投稿を表示しています。 すべての投稿を表示

2013年7月27日土曜日

RNA-seqデータの解析パイプラインを作ろう(データの視覚化)

RNA-seqのデータ(fastqファイル)をBowtieやTopHatなどでマッピングした後、bamファイルやsamファイル(場合によってはbedファイル)として出力されます。

そのデータを用いて、Cufflinksなどのリードの集計ソフトで定量値(FRKM値やRPKM値)を算出し、各転写産物の発現量を見積もると同時に、実際にゲノム上にマッピングされたリードを確認したいケースもあると思います。

そこで、今回はbam / sam / bedファイルから視覚化データを作る方法例を示したいと思います。個人でマッピングされたリードを確認する際にはIGVが便利だと思いますが、他人とデータを共有する場合では、UCSC genome browser上で可視化させたほうが都合が良いと思います。

〈スクリプトのダウンロード〉
今回は、github上で公開されているPythonスクリプトを利用します。このスクリプトはメモリの消費が少なく、複数のファイル形式(sam / bed / bowtie)に対応しているので、使い勝手が良く個人的におすすめです。

(1)rdbio-scriptsのURLに行き、右下にある「Download ZIP」をクリックし、bed2wig.pyを含む全スクリプトをダウンロードする。もしくはbed2wig.pyのURLへ移動し、スクリプトをコピペして保存する。

bed2wig.pyのスクリプト
https://github.com/daler/rdbio-scripts/blob/master/sequenceFiles/bed2wig.py
rdbio-scripts
https://github.com/daler/rdbio-scripts

〈スクリプト使用例〉
基本的な使い方は、
bed2wig.py -h
でヘルプを見ることができます。
ここでは、bamファイルをwigファイルに変換する方法について紹介したいと思います。

(1)bamファイルをソートし、samファイルに変換する。(TopHatから出力されるbamファイルは予めソートされていたと思うのでこの作業は必要ないかも。)
samtools sort accepted_hits.bam accepted_hits_sorted.bam
samtools view -F 0x0004 accepted_hits_sorted.bam > accepted_hits_sorted.sam
⇒「-F」オプションで0x0004を指定することでマッピングされなかったリードを排除しています。

(2)samファイルからwigファイルに変換する。
bed2wig.py -i accepted_hits_sorted.sam --type sam -o  accepted_hits_sorted.wig --verbose
⇒「-i」オプションでインプットファイルを指定します。「-type」オプションでインプットファイルのファイル形式を指定します。「-o」オプションでアウトプットファイルを指定します。「--verbose」オプションで進捗状況を出力させます。

(3)wigファイルの1行目を書き換える。
lessコマンドやviでファイルの中身をみてみると、1行目は「track type=wiggle_0 alwaysZero=on」となっています。このままだと使いにくいので、下記のように、適当に書き換えておきます。(もちろん、bed2wig.pyを直接書き換えても良いと思います。たぶん、269行目を書き換えればOK.。)

track name='CTRL_TopHat' description='CTRL_TopHat' type='wiggle_0' color=255,128,128 yLineMark=0.0 yLineOnOff=on visibility=2 maxHeightPixels=40:40:20 priority=1

「track name」と「description」でデータの名称を各ファイルごとに書き換えておきましょう。UCSC genome browserにアップロードする場合は、ここで指定した名称が以前アップロードしたデータと重複していると、アップロードした時に前のデータが消えてしまうので注意が必要です。

(4)ファイルを圧縮する。
UCSC genome browserではgzip (.gz), compress (.Z), or bzip2 (.bz2)の3種類の圧縮形式に対応しており、ファイルが大きくてアップロードに時間がかかりそうな場合、ファイルを圧縮しておきましょう。
bzip2 -c accepted_hits_sorted.wig > DATA.wig.bz2
⇒「-c」オプションを指定し、標準出力でファイルを出力。

〈UCSC genome browserへのアップロード〉
(1)Homeへ移動し、左側にある「Genome Browser」をクリック。
http://genome.ucsc.edu/index.html

(2)「manage custom track」ボタンをクリック。(以前にアップロードしたファイルが残っていると、「Manage Custom Tracks」ページに移動するので、その場合は「add Custom Tracks」ボタンをクリック。)

(3)「add Custom Tracks」ページに移動するので、「Paste URLs or data:」の「Or upload:」の「ファイルを選択」をクリックし、アップロードしたいwigファイルを指定する。

(4)「Submit」ボタンをクリック。ファイルがアップロードされるので、アップロードが完了するまで放置する。(アップロードしている間、ページ移動をしないこと。)

(5)アップロード完了後、「Manage Custom Tracks」ページに移動します。アップロードしたファイルは「go to Genome browser」ボタンをクリックし、Genome browser上で確認することができます。

〈Custom tracksの維持・構成の保存〉
Custom trackをアップロードしても、一週間ぐらいすると設定がリセットされて見れなくなってしまいます。そこで、Custom tracksおよびその他のTrackの構成を保存するために、UCSC genome browserが提供している「Sesson」機能を利用しましょう。

「Sesson」では、アカウントを取得することで、Custom trackの維持やTrackの構成の保存を行うことができます。また、アップロードしたCustom trackを他人が見れるように設定することもできます。

(1)Homeへ移動し、上の「My Data」→「Sessions」をクリック。
http://genome.ucsc.edu/index.html

(2)「Creat an account」をクリックし、必要なデータを入力しアカウントを作成する。

(3)Sessionにログインした後、「Save Settings」で「Save current settings as named session:」にSessionの名称を入力。

自分のアカウントにログインしていない他人にデータが見られるようにしたい場合は「allow this session to be loaded by others」のボックスにチェックをつける。

そして、「Submit」ボタンをクリックし、Sessionの保存を行う。(もちろん、Custom trackの追加やgenome browser上の構成を自分なりにアレンジして保存したい組み合わせを事前にgenome browser上で組んでおく。)

⇒Sessionを上書きしたい場合は、上書きしたいSession名を入力した後、「Submit」ボタンをクリック。上書きしてよいか確認を求められるので、「Yes」ボタンをクリック。

2013年3月20日水曜日

Integrative Genomics Viewer(IGV)を使ってみた

Integrative Genomics Viewer(IGV)とは、Standalone(自分のパソコン上)で動かせるゲノムビューワの代表格です。illuminaが提供しているGenome Studioに近いものと考えてくれればいいです。

IGVを利用する主な利点としては、

・Web上で動かす「UCSC genome browser」と比較して、動作が軽快。
⇒UCSC genome browserはサーバーを介してデータのやり取りを行うので、どうしてもタイムラグが生じる。イライラする。

・マッピング後のBAMファイルなど様々な形式のファイルを簡単に可視化できる。
⇒UCSC genome browserでは、BAMファイルの場合、サーバー上にファイルをストレージする必要があり、アップロードにも時間を要する。

・クローズドな環境でデータを見れるので、セキュリティ上安全。
⇒UCSC genome browserだとどうしてもデータを外部に提示する(アップロードする)必要があり、セキュリティ上安全とはいえないかも。

といったところでしょうか。

ではIGVのダウンロード・起動方法などについての説明に移りたいと思います。

〈IGVのダウンロード〉
(1)Integrative Genomics Viewer(IGV)のホームページへ行く。

(2)左側にある「Downloads」をクリック。

(3)「Log In」画面に移るので、
 ・初めてIGVをダウンロードする場合は、
 (ⅰ)「To use IGV, registration is required. Click here to register」で「Click here」をクリック。

 (ⅱ)Registrationに移るので、「Name」「Email」「Organization」を記入後、「Agree」をクリック。

 ・Registrationを済ませている場合は、
 (ⅰ)「email address」を記入後、「Login」ボタンをクリック。

(3)「Downloads」画面に移るので、「Binary Distribution」のIGV_2.2.11.zipをクリックし、IGV本体をダウンロードする。

(4)また「igvtools」のigvtools_2.2.2.zip     includes .genome files (225 MB) をクリックし、Genome情報をあわせてダウンロードしておく。

〈IGVの起動方法〉
(1)Zipファイルを解凍した中に、「igv.bat」というファイルがあり、そのファイルをクリックするとIGVが起動します。

〈使用するゲノムを導入〉
デフォルトではヒトゲノムである「hg18」しか入っていないので、これ以外のゲノムを使いたい場合は自身でIGVの中にファイルをインポートする必要があります。

そこで、先ほどダウンロードしてきた「igvtools」中のファイルを使います。
「igvtools」にはさまざまなゲノム情報のファイルが入っており、このファイルをIGVの中にインポートすることで目的のゲノムを使用することができるようになります。

(1)「igvtools」を解凍する。

(2)「igvtools」→「genomes」内にさまざまな生物種のゲノム情報が入っています。(場所を確認)

(3)IGVを起動し、メニューバーの「Genomes」→「Load genome from File…」をクリック。

(4)使用したいゲノム情報をクリック(hg19.genome等のファイル)

(5)「Loading Genome…」というウインドウが出るのでしばらく待つ。

(6)ロードが終わった後、メニューバーの右端下のプルダウンメニューから、先ほどインポートしたゲノムが選択できるようになっている。

〈BAMファイルの可視化〉
IGVでBAMファイルを可視化するにはいくつか条件があります。

・「Chromosome」および「Start Position」でソートする必要あり。
・マッピングしたファイル(.bam)およびIndexファイル(bam.bai)が必要。

「Tophat」でマッピングした場合、BAMファイルはソートされた状態で生成されますが、「Bowtie」でマッピングした場合、生成されるSAMファイルはソートされていないので注意が必要です。

SAMtoolsを使って、sortとindexの生成を行います。
コマンドは以下のとおり。(TophatからBAMファイルを生成した場合は(3)から。)

(1)SAMファイルからBAMファイルへの変換
$samtools view -bS INPUT_file.sam > OUTPUT_file.bam
⇒標準出力により、BAMファイルは生成されます。

(2)BAMファイルのソート
$ samtools sort INPUT_file.bam OUTPUT_file
⇒アウトプットのファイルには自動的に「.bam」の拡張子がつくのでコマンド内で記述する必要はありません。

(3)BAMファイルのindex生成
$ samtools index INPUT_file.bam(sort済み)
⇒今度はOUTPUTファイルの名前を指定する必要はありません。上記の例では、「INPUT_file.bam.bai」というindexファイルが作られます。

〈IGVでBAMファイルを可視化〉
(1)IGVを起動。(「igv.bat」ファイルをクリック)

(2)メニューから「File」→「Load from File…」をクリック。

(3)見たいBAMファイル(.bam)をクリック。
⇒このとき、あわせてindexファイル(.bam.bai)を開く必要はない。BAMファイルと同じ場所にindexファイルがあれば問題なくBAMファイルを開くことができる。

2012年12月24日月曜日

RNA-seqのデータをUCSC genome browser上で視覚化する(3)

前回作成したcustom trackや他の既存のトラックの組み合わせを保存しておきたいといったとき、どうすればよいか。そんなとき、UCSC genome browserが提供する「Sessions」と呼ばれるトラック管理機能を使うと便利です。

〈主なステップ〉
1. Track linesの記述
2. データのアップロード
3. データの管理方法

3. データの管理方法
(1)UCSC genome browserのトップページから、上にあるメニューの「Session」をクリック

(2)Sign in to UCSC Genome Bioinformaticsの「Login」
⇒アカウントを持っていない場合は、「Creat an account」からアカウントを作成してください。

(3)Session ManagementのSave Settingsにある「Save current settings as named session」でセッションの名前を入力し、「Submit」をクリック
⇒他の人と、そのセッションを共有したい場合は「allow this session to be loaded by others」にチェックを入れると、他の人もそのセッションを閲覧可能になる。

(4)My sessionsに先ほど作ったセッションが登録される。link to sessionの「Browser」をクリックすることで保存したトラックを閲覧でき、share with others?にチェックを入れることで他の人の閲覧が可能となる。
⇒各セッションのBrowserのURLから直接自分の組み合わせたトラックにアクセスできるので、そのURLを他の人に教えることで、他の人とデータを共有することができる。(セキュリティ上の問題などがあるので、注意して使用する必要あり。)

ちなみに、Custom trackとしてアップロードしたデータは、4ヶ月以上アクセスしないで放置していると、データが消えてしまうので注意してください。

〈参考〉
■UCSC genome browser -Sessionの使用方法について
http://genome.ucsc.edu/goldenPath/help/hgSessionHelp.html


2012年12月23日日曜日

RNA-seqのデータをUCSC genome browser上で視覚化する(2)

UCSC genome browser上で可視化できるデータのフォーマットとして、bedGraph, GTF, BED, WIG, bigwig, BAMが主なものとして挙げられます。
今回は、前回作成した「WIG」と呼ばれる形式のファイルを利用してRNA-seqのデータを可視化してみたいと思います。

〈主なステップ〉
1. Track linesの記述
2. データのアップロード
3. データの管理方法

1. Track linesの記述
まず注意しなければいけない点は、Tophatなどのマッピングソフトから得られた「BAMファイル」やGENCODEなどからダウンロードしてきた「GTFファイル」などをアップロードする前に、ファイルの冒頭(最初の一行目)に「Track lines」と呼ばれる一種の但し書きのようなものを書く必要があるということです。

これは、UCSC genome browser上にアップロードしたデータのファイル形式(BAM, GTFなど)、Genome browser上で表示させる時のTrackの名前・色などの基本情報を付与するためにあります。

では具体的にどのように記述すればよいか見て行きましょう。
まずは、Track linesの定義についてよく利用するものを中心に説明していきたいと思います。

〈Track lines〉
・track name="track_label"
 Genome browserのウインドウの右側に位置するラベル名を定義。
・description="center_label"
 Genome browserのウインドウの中央に位置するラベル名を定義。
 60文字以内という字数制限がある。
 空白はできるだけ使わずに「ハイフン"_"」を使ったほうが良い。
・type="track_type"
 アップロードするファイルの形式を定義。
 WIGファイルなら「Wiggle_0」、BEDファイルなら「bed」と記載。 
・visibility=number
 annotation trackのデフォルトの表示モードを定義。
 0-4の数字に表示モードが対応しており、
 「0 -hide」「1-dense」「2-full」「3-pack」「4-squish」となっている。
 定義しない場合では、「1-dense」が自動的に選択される。
・color=RRR,GGG,BBB
 annotation trackの色を定義。
 「コンマ","」で区切られた0-255の幅を持つRGB valuesによって色が定義されている。
 定義しない場合では、「0,0,0」(黒)が自動的に選択される。
 色の組み合わせを考えるときの参考として、
 MUDCUBE -COLOR SPHEREなどのサイトを活用するといいかも。
 (16進法による色の記述になっているので、RGB変換する必要あり。)
・colorByStrand=RRR,GGG,BBB,RRR,GGG,BBB
 Genome上の「+鎖」「-鎖」を区別して色分けできる。
 
2. データのアップロード
(1)UCSC genome browserのトップページから上のメニューの「Genome」、もしくは左側のメニューの「Genome Browser」をクリック
(2)「manage custom tracks」をクリック
(3)「add custom tracks」をクリック
(4)Paste URLs or data:にある「ファイルを選択」をクリックし、
アップロードしたいファイルを選択
(5)「Submit」をクリックし、ファイルをアップロード(しばらく時間がかかるので放置)
(6)アップロード終了後、「go to genome browser」をクリックすると、
アップロードしたcustom trackを確認できます。


〈具体例〉
NONCODE v3.0のlncRNA_humanの「BEDファイル」のデータをUCSC genome browser上で可視化してみよう。
■使用したデータ
http://www.noncode.org/NONCODERv3/datadownload/lncRNA_human.zip

例1:
まずは、最低限の情報として「track name」「description」「type」のみを1行目に記述したファイルを作成してみる。
track name="NONCODE_lncRNA_human" description="NONCODE" type="bed"
chr7 89010556 89010766 n123 1000 +
chr7 18847664 18847902 n125 1000 -
chr7 92600316 92600610 n127 1000 -
chr12 6619387 6619717 n1315 1000 +
…

例2:
例1と少し内容を追加して、アップロードしたトラックの違いを比較してみる。下記のように「visibility」「color」の記述を追加すると、デフォルトの表示モードがはじめから「Full」の状態になり、トラックに色がついていることが確認できる。
track name="NONCODE_lncRNA_human" description="NONCODE" type="bed" visibility=2 color=0,142,247
chr7 89010556 89010766 n123 1000 +
chr7 18847664 18847902 n125 1000 -
chr7 92600316 92600610 n127 1000 -
chr12 6619387 6619717 n1315 1000 +
…


〈参考〉
■UCSC genome browserへのアップロードに関する注意事項
http://genome.ucsc.edu/goldenPath/help/customTrack.html
■MUDCUBE -COLOR SPHERE
http://mudcu.be/sphere/
■RGB変換
http://www.kitaq.net/lib/rgb/



2012年12月17日月曜日

RNA-seqのデータをUCSC genome browser上で視覚化する(1)

Bowtieなどのマッピングソフトから出力される「bamファイル」を「Wigファイル」に変換し、UCSC genome browserにアップロードできるファイル形式にする。

(1)bamToBed
ダウンロード先:http://code.google.com/p/bedtools/downloads/list

〈コマンド〉
$ bamToBed -i  Input_name > Output_name
$ sortBed -i Input_name > Output_name

(2)bed2Wig
ダウンロード先:不明

〈コマンド〉
$ cut -f 1-3 Input_name > Output_name
$ perl bed2Wig.pl Input_name > Output_name

〈スクリプト〉
ファイル名:bed2Wig.pl

内容:
#!/usr/local/bin/perl
use strict;
my $m;
while (<>) {
    my @tok = split "\t";
    for ($tok[1]..$tok[2]) {
$m->{$tok[0]}{$_}++
    }
}
my $track = shift || "name";
my $desc = shift || "desc";
#track name / description: 任意の名前に入力, color: 任意の色を指定
print "track name=\"$track\" description=\"$desc\" type=\"wiggle_0\" color=250,180,10 yLineMark=0.0 yLineOnOff=on visibility=2 maxHeightPixels=40:40:20 priority=1\n";
for my $chr (sort keys %$m) {
    my $prev = 0;
    for my $pos (sort {$a<=>$b} keys %{$m->{$chr}}) {
    if ($pos - $prev > 1) {
   print "fixedStep chrom=$chr start=$pos step=1\n";
    }
    print $m->{$chr}{$pos}, "\n";
    $prev = $pos;
    }
}