そこで、Linuxの改行コードに変換するcommandを紹介します。
dos2unix ファイル名でLinuxの改行コードに変換可能です。データは上書きされます。
2013年9月14日土曜日
WindowsとLinux間での改行コードの違いによる問題
改行コードの違いから、Perlなどのスクリプトで処理した時、予期せぬところに改行が入ったりすることがあります。
そこで、Linuxの改行コードに変換するcommandを紹介します。
そこで、Linuxの改行コードに変換するcommandを紹介します。
Illumina Adapter Sequences Letter: アダプター配列等の情報を探す（illumina）
fastqファイルのpreprocessingを段階で、アダプター配列を除去したい場合があると思います（small RNAのRNA-seqやCLIP-seqなど3'末端にリンカー配列を結合させたケースなど）。
しかし残念ながら、論文のMaterial&Methodsをみてもそのアダプターの配列が記載されていないことがあります。
アダプター配列がわからない場合、FastQCなどでリード中に高頻度で出現する配列をチェックして、それらの配列を除去するというのも１つの手ですが、単純にメーカーのサイトから配列情報を入手することも可能です。
Illuminaの場合、下記のサイトから「Illumina Adapter Sequences Letter」と呼ばれるアダプター配列等の配列情報が記載されたファイルをダウンロードできます。
http://support.illumina.com/downloads/illumina_adapter_sequences_letter.ilmn
このように、論文中に配列情報が記載されていなくてもサンプルプレップで使用したキットがわかっていればアダプター配列を特定することが可能です。
しかし残念ながら、論文のMaterial&Methodsをみてもそのアダプターの配列が記載されていないことがあります。
アダプター配列がわからない場合、FastQCなどでリード中に高頻度で出現する配列をチェックして、それらの配列を除去するというのも１つの手ですが、単純にメーカーのサイトから配列情報を入手することも可能です。
Illuminaの場合、下記のサイトから「Illumina Adapter Sequences Letter」と呼ばれるアダプター配列等の配列情報が記載されたファイルをダウンロードできます。
http://support.illumina.com/downloads/illumina_adapter_sequences_letter.ilmn
このように、論文中に配列情報が記載されていなくてもサンプルプレップで使用したキットがわかっていればアダプター配列を特定することが可能です。
liftOver: ゲノム座標を異なるバージョンのゲノム座標に対応付けする
hg18のゲノム座標のデータを持っていて、hg19のゲノム座標に変換したい。そんな場面があるかと思います。
こういった異なるバージョン間でのゲノム座標の対応付けは、UCSCが提供している「liftOver」が便利です。
〈liftOverのダウンロード〉
（１）下記のURLからliftOverをダウンロード
・32bit版（Linux）
http://hgwdev.cse.ucsc.edu/~kent/exe/linux/liftOver.gz
・64bit版（Linux）
http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver
〈liftOverのインストール〉
（１）ダウンロードしたliftOverのファイル（32bit版では解凍する必要あり）に実行権限を与える
・変換テーブル(例えば、「hg18⇒hg19」では、hg18ToHg19.over.chain)
（１）UCSC genome browserの「Downloads」から、特定の種（e.g. human）をクリックする。
例えば、http://hgdownload.soe.ucsc.edu/downloads.html#human
（２）移動先から、Inputファイルのゲノム座標のバージョンの場所まで移動
⇒下を見ていくと、例えば、「Mar. 2006 (hg18)」など各バージョンのデータが箇条書きになって並んでいるはず。
（３）「LiftOver files」をクリック。
（４）特定の変換テーブルを選択、クリックしデータをダウンロード。
・Bedファイル（変換したいゲノム座標の情報）
⇒各自で用意。
〈liftOver実行例〉
・１つ目: 座標軸を変換したいインプットファイル（Bedファイル）の名前・ディレクトリを指定
・２つ目: 変換テーブルの名前・ディレクトリの指定
・３つ目: 変換後のアウトプットファイルの名前・ディレクトリの指定
・４つ目: 変換できなかったものの出力先
〈参考リンク〉
・NGS Surfer's Wiki - ゲノムのバージョンの違う情報を対応付ける
http://cell-innovation.nig.ac.jp/wiki/tiki-index.php?page=%E3%82%B2%E3%83%8E%E3%83%A0%E3%81%AE%E3%83%90%E3%83%BC%E3%82%B8%E3%83%A7%E3%83%B3%E3%81%AE%E9%81%95%E3%81%86%E6%83%85%E5%A0%B1%E3%82%92%E5%AF%BE%E5%BF%9C%E4%BB%98%E3%81%91%E3%82%8B
こういった異なるバージョン間でのゲノム座標の対応付けは、UCSCが提供している「liftOver」が便利です。
〈liftOverのダウンロード〉
（１）下記のURLからliftOverをダウンロード
・32bit版（Linux）
http://hgwdev.cse.ucsc.edu/~kent/exe/linux/liftOver.gz
・64bit版（Linux）
http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver
〈liftOverのインストール〉
（１）ダウンロードしたliftOverのファイル（32bit版では解凍する必要あり）に実行権限を与える
sudo chmod +x liftOver（２）パスの通っているディレクトリ（e.g. /bin）に移動
sudo mv liftOver /bin〈必要となるファイルの準備〉
・変換テーブル(例えば、「hg18⇒hg19」では、hg18ToHg19.over.chain)
（１）UCSC genome browserの「Downloads」から、特定の種（e.g. human）をクリックする。
例えば、http://hgdownload.soe.ucsc.edu/downloads.html#human
（２）移動先から、Inputファイルのゲノム座標のバージョンの場所まで移動
⇒下を見ていくと、例えば、「Mar. 2006 (hg18)」など各バージョンのデータが箇条書きになって並んでいるはず。
（３）「LiftOver files」をクリック。
（４）特定の変換テーブルを選択、クリックしデータをダウンロード。
・Bedファイル（変換したいゲノム座標の情報）
⇒各自で用意。
〈liftOver実行例〉
liftOver INPUT.bed hg18ToHg19.over.chain OUTPUT.bed ERROR.bedコマンドライン引数が４つある。
・１つ目: 座標軸を変換したいインプットファイル（Bedファイル）の名前・ディレクトリを指定
・２つ目: 変換テーブルの名前・ディレクトリの指定
・３つ目: 変換後のアウトプットファイルの名前・ディレクトリの指定
・４つ目: 変換できなかったものの出力先
〈参考リンク〉
・NGS Surfer's Wiki - ゲノムのバージョンの違う情報を対応付ける
http://cell-innovation.nig.ac.jp/wiki/tiki-index.php?page=%E3%82%B2%E3%83%8E%E3%83%A0%E3%81%AE%E3%83%90%E3%83%BC%E3%82%B8%E3%83%A7%E3%83%B3%E3%81%AE%E9%81%95%E3%81%86%E6%83%85%E5%A0%B1%E3%82%92%E5%AF%BE%E5%BF%9C%E4%BB%98%E3%81%91%E3%82%8B
2013年8月10日土曜日
grep: and / or条件での検索
Linuxのコマンドの１つである「grep」でAND条件・OR条件での検索に関するメモ。
（１）AND検索
パイプ「|」でつなぐことでgrepコマンドのAND検索が可能になります。
grepのあとの"ダブルクオーテーション"内に検索したいキーワード（もしくは正規表現）を記述しますが、キーワード（例えばA\|B）の間を「\|」(バックスラッシュとバーティカルバー)で区切ることで、OR検索が可能になります（例えばAもしくはB）。
（１）AND検索
パイプ「|」でつなぐことでgrepコマンドのAND検索が可能になります。
grep "chr6:125111" INPUT.txt | grep "1:0:0" > OUTPUT.txt（２）OR検索
grepのあとの"ダブルクオーテーション"内に検索したいキーワード（もしくは正規表現）を記述しますが、キーワード（例えばA\|B）の間を「\|」(バックスラッシュとバーティカルバー)で区切ることで、OR検索が可能になります（例えばAもしくはB）。
grep "chr6:125111\|chr6:125112" INPUT.txt > OUTPUT.txt
2013年8月3日土曜日
sortBed: bedファイルをソートする
bedファイルをsortコマンドでソートしようとしたことがあるのですがこれが意外と難しい(´・ω・`)
そこで便利なのが、BedtoolsのsortBedコマンドです。
〈参考文献〉
The BEDTools suite
http://bedtools.readthedocs.org/en/latest/content/tools/sort.html
そこで便利なのが、BedtoolsのsortBedコマンドです。
sortBed -i /path/to/INPUT.bed > OUTPUT.bedBED, GFF, VCFの３種類のファイルに対応しており、これらのファイルをソートすることができます。
〈参考文献〉
The BEDTools suite
http://bedtools.readthedocs.org/en/latest/content/tools/sort.html
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）に対応しているので、使い勝手が良く個人的におすすめです。
（１）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
〈スクリプト使用例〉
基本的な使い方は、
ここでは、bamファイルをwigファイルに変換する方法について紹介したいと思います。
（１）bamファイルをソートし、samファイルに変換する。（TopHatから出力されるbamファイルは予めソートされていたと思うのでこの作業は必要ないかも。）
（２）samファイルからwigファイルに変換する。
（３）wigファイルの１行目を書き換える。
lessコマンドやviでファイルの中身をみてみると、１行目は「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にアップロードする場合は、ここで指定した名称が以前アップロードしたデータと重複していると、アップロードした時に前のデータが消えてしまうので注意が必要です。
（４）ファイルを圧縮する。
UCSC genome browserではgzip (.gz), compress (.Z), or bzip2 (.bz2)の３種類の圧縮形式に対応しており、ファイルが大きくてアップロードに時間がかかりそうな場合、ファイルを圧縮しておきましょう。
〈UCSC genome browserへのアップロード〉
（１）Homeへ移動し、左側にある「Genome Browser」をクリック。
http://genome.ucsc.edu/index.html
（２）「manage custom track」ボタンをクリック。（以前にアップロードしたファイルが残っていると、「Manage Custom Tracks」ページに移動するので、その場合は「add Custom Tracks」ボタンをクリック。）
（３）「add Custom Tracks」ページに移動するので、「Paste URLs or data:」の「Or upload:」の「ファイルを選択」をクリックし、アップロードしたいwigファイルを指定する。
（４）「Submit」ボタンをクリック。ファイルがアップロードされるので、アップロードが完了するまで放置する。（アップロードしている間、ページ移動をしないこと。）
（５）アップロード完了後、「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を他人が見れるように設定することもできます。
（１）Homeへ移動し、上の「My Data」→「Sessions」をクリック。
http://genome.ucsc.edu/index.html
（２）「Creat an account」をクリックし、必要なデータを入力しアカウントを作成する。
（３）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」ボタンをクリック。
そのデータを用いて、Cufflinksなどのリードの集計ソフトで定量値（FRKM値やRPKM値）を算出し、各転写産物の発現量を見積もると同時に、実際にゲノム上にマッピングされたリードを確認したいケースもあると思います。
そこで、今回はbam / sam / bedファイルから視覚化データを作る方法例を示したいと思います。個人でマッピングされたリードを確認する際にはIGVが便利だと思いますが、他人とデータを共有する場合では、UCSC genome browser上で可視化させたほうが都合が良いと思います。
〈スクリプトのダウンロード〉
今回は、github上で公開されているPythonスクリプトを利用します。このスクリプトはメモリの消費が少なく、複数のファイル形式（sam / bed / bowtie）に対応しているので、使い勝手が良く個人的におすすめです。
（１）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ファイルに変換する方法について紹介したいと思います。
（１）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を指定することでマッピングされなかったリードを排除しています。
（２）samファイルからwigファイルに変換する。
bed2wig.py -i accepted_hits_sorted.sam --type sam -o accepted_hits_sorted.wig --verbose⇒「-i」オプションでインプットファイルを指定します。「-type」オプションでインプットファイルのファイル形式を指定します。「-o」オプションでアウトプットファイルを指定します。「--verbose」オプションで進捗状況を出力させます。
（３）wigファイルの１行目を書き換える。
lessコマンドやviでファイルの中身をみてみると、１行目は「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にアップロードする場合は、ここで指定した名称が以前アップロードしたデータと重複していると、アップロードした時に前のデータが消えてしまうので注意が必要です。
（４）ファイルを圧縮する。
UCSC genome browserではgzip (.gz), compress (.Z), or bzip2 (.bz2)の３種類の圧縮形式に対応しており、ファイルが大きくてアップロードに時間がかかりそうな場合、ファイルを圧縮しておきましょう。
bzip2 -c accepted_hits_sorted.wig > DATA.wig.bz2⇒「-c」オプションを指定し、標準出力でファイルを出力。
〈UCSC genome browserへのアップロード〉
（１）Homeへ移動し、左側にある「Genome Browser」をクリック。
http://genome.ucsc.edu/index.html
（２）「manage custom track」ボタンをクリック。（以前にアップロードしたファイルが残っていると、「Manage Custom Tracks」ページに移動するので、その場合は「add Custom Tracks」ボタンをクリック。）
（３）「add Custom Tracks」ページに移動するので、「Paste URLs or data:」の「Or upload:」の「ファイルを選択」をクリックし、アップロードしたいwigファイルを指定する。
（４）「Submit」ボタンをクリック。ファイルがアップロードされるので、アップロードが完了するまで放置する。（アップロードしている間、ページ移動をしないこと。）
（５）アップロード完了後、「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を他人が見れるように設定することもできます。
（１）Homeへ移動し、上の「My Data」→「Sessions」をクリック。
http://genome.ucsc.edu/index.html
（２）「Creat an account」をクリックし、必要なデータを入力しアカウントを作成する。
（３）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年7月21日日曜日
RNA-seqデータの解析パイプラインを作ろう（TopHat-Cuffilnks）
RNA-seqデータの解析例をメモ。
あくまで一例なので、ここで取り上げた方法がすべてのデータに対して良いとは限りません。
一つ一つコマンドを入力していってもいいのですが、手間なのでシェルスクリプトを活用し複数の処理（コマンド）をまとめて行った例を示します（いわゆる、バッチ処理）。
以下では、RNA-seqデータの解析例をシェルスクリプトを示しました。
・Illumina CASAVA-1.8 FASTQ Filter
・FASTX-Toolkit 0.0.13
・Bowtie 0.12.7
・SAM tools 0.1.18
・TopHat 2.0.4
・Cufflinks 2.1.1
〈使用したファイル〉
1) Ribosomal_RNA
Ribosomal RNA（18S, 28S）のIndexファイル
⇒NCBIなどからFasta形式で配列情報を取得し、そのファイルからBowtie-buildコマンドでIndexファイルを事前に用意しておき、Bowtieのディレクトリの中の「indexes」ディレクトリにダウンロードしたIndexファイルを入れておく。
2) hg19
Human genome（hg19）のIndexファイル
http://bowtie-bio.sourceforge.net/index.shtml
⇒Bowtieのサイトの右側に、「Pre-build indexes」の項目があり、様々な種のゲノム配列情報を取得することができる。Fastqファイルをゲノムにマッピング際に必要なので、事前に該当のデータをダウンロードしておき、Bowtieのディレクトリの中の「indexes」ディレクトリにダウンロードしたIndexファイルを入れておく。
3) Refseq_genes_May_14_2013.gtf
Illumina igenomesの「transcript.gtf」ファイル
http://cufflinks.cbcb.umd.edu/igenomes.html
（１）一行目に「#!/bin/sh」を書く
（２）スクリプトのファイルに実行権限を与える
の２点が必要になります。
（２）に関しては、以下のcommandを実行し、事前に実行権限を付与しておきましょう。
（１）クオリティチェック～マッピング(ファイル名: RNA-seq_pipeline.sh)
順を追って説明いきます。
1行目の「#!/bin/sh」はこのファイルがシェルスクリプトであることを明示するためのものであり、指定したパス「/bin/sh」のプログラムによって該当のファイルを実行させるのに必要です。
2行目はループ関数である「for」を用いて、複数のFastqファイルに対して、doからdoneの間の処理を繰り返し行わせています。ここで、「$@」は1つ以上のコマンドライン引数（SRRxxxxxxのようなファイル名）を受け取り、「file」に代入しています。
5行目では、「Illumina CASAVA-1.8 FASTQ Filter」のソフトを利用し、CASAVA filterをパスしなかったクオリティの低いリードを除去します（Illumina以外のプラットフォームから得られたFastqファイルやNCBI GEOやDDBJ SRAからダウンロードしてきたFastqファイルに対しては不要）。
HiseqなどのIlluminaのシークエンサーから得られたFastqファイルには、CASAVA filterをパスしたかどうかの目印が各リードに与えられていますので（フィルターをパスしていれば「N」、フィルターをパスしていなければ「Y」）、ソフトを利用しフィルターをパスしたリードのみ抽出します。
「-N」オプションによりフィルターをパスしたリードのみを出力させ、「-o」オプションにより、指定したディレクトリに出力したFastqファイルを保存させます。
10行目では、さらにクオリティの低いリードのフィルタリング、及びクオリティの低い塩基のトリミングを行なっています。具体的には、Phred quality score 20未満が80%以上のリードを除去、Phred quality score 20未満の末端を除去、配列長が10bp未満のリードを除去を行なっています。
ここで、Phred quality scoreとは、Fastqファイル中に記載されているリードのクオリティを示す指標のことです。Phred quality scoreは以下の式によって表されます。
Phred quality score (Q) = -10 log10P
Pは「Possibility」、つまりはシークエンサーから読まれた各塩基の誤差率を示しています。例えば、P = 20であれば、誤差率は1%であり、100回に1回読み間違いが起こることを意味しています。そのため、リードのクオリティチェックでは、リードの各塩基のクオリティが「P >= 20」もしくは「P >= 30」ぐらいになるようにフィルターをかけるとよいのではないかと思います（シークエンサーで読んだリード長に応じてクライテリアを設定）。
13行目では、Ribosomal RNA由来のリードを排除します。「--un」オプションは指定したディレクトリにマッピングされなかったリードをFastqファイルで出力させます。また、Bowtieでは、コマンドライン引数は「2つ」指定する必要があり、１つ目は、Fastqファイル中のリードをマッピングするリファレンス配列（ゲノム配列・その他の配列）を指定し、2つ目は、マッピングしたいFastqファイルを指定します。上記では、コマンドライン引数の1つ目が「Ribosomal_RNA」、2つ目が「"$file"_filtered.fastq」に該当します。また、標準出力としてマッピングされたリード（Ribosomal RNA由来のリード）を「/dev/null」（ゴミ箱のようなもの）へリダイレクトするように指定しています。
「-p」オプションでは、パソコンのコア数を指定しています。デフォルトでは「-p 2」となっており、2コア分で計算を行なっています。CPUのコア数が2よりも多い場合は、自分のパソコンのコア数に応じて指定したほうがよいです。その分、計算速度も向上します。
16行目で、いよいよリファレンス配列へのマッピングを行なっていきます。対象がヒトである場合、パラメータはいじくらなくてもヒト用に最適化されている（はず）なので、デフォルトでとりあえずマッピングしてます。
「-o」オプションで出力先のディレクトリを指定します（ディレクトリがない場合は生成してくれます）。
「-G」オプションで使いたいアノテーション情報（gtfファイル）を指定します。また、コマンドライン引数についてはBowtieの時と同様です。
（２）FRPM値(RPKM値)の算出
Cuffilnksのコマンドである「Cuffdiff」を用いてFPKM値（RPKM値）の算出と発現変動遺伝子の同定を行います。
コマンドライン引数は3つ指定する必要があります。1つ目は、アノテーション情報（gtfファイル）を指定し、2つ目と3つ目には比較したいデータ（bamファイル）をそれぞれ指定します。Duplicateのデータ（n >= 2）がある場合は、各データ（bamファイル）を「,」カンマで区切り指定することができます（比較するデータをn = 1でしか取っていない場合、有意差を計算することができません。そのため、この場合Fold-change（例えば2以上）をもとに発現変動遺伝子を絞り込む必要があります）。
〈シェルスクリプト実行例〉
/path/to/RNA-seq_pipeline.sh SRRxxxxxx SRRyyyyyy SRRzzzzzz
⇒シェルスクリプトを呼び出し、コマンドライン引数として処理したいFastqファイルの名前を指定。
あくまで一例なので、ここで取り上げた方法がすべてのデータに対して良いとは限りません。
一つ一つコマンドを入力していってもいいのですが、手間なのでシェルスクリプトを活用し複数の処理（コマンド）をまとめて行った例を示します（いわゆる、バッチ処理）。
以下では、RNA-seqデータの解析例をシェルスクリプトを示しました。
#!/bin/sh for file in "$@" do #Illumina CASAVA Filter [Y]を除去 fastq_illumina_filter -N -o /path/to/"$file"_casava.fastq /path/to/"$file".fastq #Phred quality score 20未満が80%以上のリードを除去 #Phred quality score 20未満の末端を除去 #配列長が10bp未満のリードを除去 fastq_quality_filter -Q33 -q 20 -p 80 -i /path/to/"$file"_casava.fastq | fastq_quality_trimmer -Q33 -t 20 -l 10 -o /path/to/"$file"_filtered.fastq #Ribosomal RNA（rRNA）由来のリードの除去 bowtie -p 8 --un /path/to/"$file"_norrna.fastq /path/to/bowtie-X.XX.X/indexes/Ribosomal_RNA /path/to/"$file"_filtered.fastq > /dev/null #Reference genome / Transcriptomeへのマッピング tophat -p 8 -o /path/to/"$file"_out -G /path/to/Refseq_genes_May_14_2013.gtf /path/to/bowtie-X.XX.X/indexes/hg19 /path/to/"$file"_norrna.fastq done
cuffdiff -p 8 -o /path/to/cuffdiff_out /path/to/Refseq_genes_May_14_2013.gtf <CTRL_rep1.bam,CTRL_rep2.bam…> <KD_rep1.bam,KD_rep2.bam…>〈使用したソフトウェア〉
・Illumina CASAVA-1.8 FASTQ Filter
・FASTX-Toolkit 0.0.13
・Bowtie 0.12.7
・SAM tools 0.1.18
・TopHat 2.0.4
・Cufflinks 2.1.1
〈使用したファイル〉
1) Ribosomal_RNA
Ribosomal RNA（18S, 28S）のIndexファイル
⇒NCBIなどからFasta形式で配列情報を取得し、そのファイルからBowtie-buildコマンドでIndexファイルを事前に用意しておき、Bowtieのディレクトリの中の「indexes」ディレクトリにダウンロードしたIndexファイルを入れておく。
2) hg19
Human genome（hg19）のIndexファイル
http://bowtie-bio.sourceforge.net/index.shtml
⇒Bowtieのサイトの右側に、「Pre-build indexes」の項目があり、様々な種のゲノム配列情報を取得することができる。Fastqファイルをゲノムにマッピング際に必要なので、事前に該当のデータをダウンロードしておき、Bowtieのディレクトリの中の「indexes」ディレクトリにダウンロードしたIndexファイルを入れておく。
3) Refseq_genes_May_14_2013.gtf
Illumina igenomesの「transcript.gtf」ファイル
http://cufflinks.cbcb.umd.edu/igenomes.html
⇒上記のURLから「Homo sapiens - UCSC - hg19」を選択、該当のファイル（20GBを超える）をダウンロード。その中に「transcript.gtf」があり、このファイルの中にRefseqのアノテーション情報が記載されている。
⇒目的に応じて、gencode / Ensemblなどの他のgtfファイルを利用することも可。ただ、多くのIsoformが登録されているので、どこまで正確に各Isoformごとの定量値を見積もることができるのか疑問。（そもそも、Isoformごとの正確な定量値を必要なのかどうかにもよる。）
〈シェルスクリプトを利用する際の注意点〉
rootユーザー以外が上記のシェルスクリプトを使用する場合、（１）一行目に「#!/bin/sh」を書く
（２）スクリプトのファイルに実行権限を与える
の２点が必要になります。
（２）に関しては、以下のcommandを実行し、事前に実行権限を付与しておきましょう。
$ chmod +x /path/to/RNA-seq_pipeline.sh〈シェルスクリプトの内容について〉
（１）クオリティチェック～マッピング(ファイル名: RNA-seq_pipeline.sh)
順を追って説明いきます。
1行目の「#!/bin/sh」はこのファイルがシェルスクリプトであることを明示するためのものであり、指定したパス「/bin/sh」のプログラムによって該当のファイルを実行させるのに必要です。
2行目はループ関数である「for」を用いて、複数のFastqファイルに対して、doからdoneの間の処理を繰り返し行わせています。ここで、「$@」は1つ以上のコマンドライン引数（SRRxxxxxxのようなファイル名）を受け取り、「file」に代入しています。
5行目では、「Illumina CASAVA-1.8 FASTQ Filter」のソフトを利用し、CASAVA filterをパスしなかったクオリティの低いリードを除去します（Illumina以外のプラットフォームから得られたFastqファイルやNCBI GEOやDDBJ SRAからダウンロードしてきたFastqファイルに対しては不要）。
HiseqなどのIlluminaのシークエンサーから得られたFastqファイルには、CASAVA filterをパスしたかどうかの目印が各リードに与えられていますので（フィルターをパスしていれば「N」、フィルターをパスしていなければ「Y」）、ソフトを利用しフィルターをパスしたリードのみ抽出します。
「-N」オプションによりフィルターをパスしたリードのみを出力させ、「-o」オプションにより、指定したディレクトリに出力したFastqファイルを保存させます。
10行目では、さらにクオリティの低いリードのフィルタリング、及びクオリティの低い塩基のトリミングを行なっています。具体的には、Phred quality score 20未満が80%以上のリードを除去、Phred quality score 20未満の末端を除去、配列長が10bp未満のリードを除去を行なっています。
ここで、Phred quality scoreとは、Fastqファイル中に記載されているリードのクオリティを示す指標のことです。Phred quality scoreは以下の式によって表されます。
Phred quality score (Q) = -10 log10P
Pは「Possibility」、つまりはシークエンサーから読まれた各塩基の誤差率を示しています。例えば、P = 20であれば、誤差率は1%であり、100回に1回読み間違いが起こることを意味しています。そのため、リードのクオリティチェックでは、リードの各塩基のクオリティが「P >= 20」もしくは「P >= 30」ぐらいになるようにフィルターをかけるとよいのではないかと思います（シークエンサーで読んだリード長に応じてクライテリアを設定）。
13行目では、Ribosomal RNA由来のリードを排除します。「--un」オプションは指定したディレクトリにマッピングされなかったリードをFastqファイルで出力させます。また、Bowtieでは、コマンドライン引数は「2つ」指定する必要があり、１つ目は、Fastqファイル中のリードをマッピングするリファレンス配列（ゲノム配列・その他の配列）を指定し、2つ目は、マッピングしたいFastqファイルを指定します。上記では、コマンドライン引数の1つ目が「Ribosomal_RNA」、2つ目が「"$file"_filtered.fastq」に該当します。また、標準出力としてマッピングされたリード（Ribosomal RNA由来のリード）を「/dev/null」（ゴミ箱のようなもの）へリダイレクトするように指定しています。
「-p」オプションでは、パソコンのコア数を指定しています。デフォルトでは「-p 2」となっており、2コア分で計算を行なっています。CPUのコア数が2よりも多い場合は、自分のパソコンのコア数に応じて指定したほうがよいです。その分、計算速度も向上します。
16行目で、いよいよリファレンス配列へのマッピングを行なっていきます。対象がヒトである場合、パラメータはいじくらなくてもヒト用に最適化されている（はず）なので、デフォルトでとりあえずマッピングしてます。
「-o」オプションで出力先のディレクトリを指定します（ディレクトリがない場合は生成してくれます）。
「-G」オプションで使いたいアノテーション情報（gtfファイル）を指定します。また、コマンドライン引数についてはBowtieの時と同様です。
（２）FRPM値(RPKM値)の算出
Cuffilnksのコマンドである「Cuffdiff」を用いてFPKM値（RPKM値）の算出と発現変動遺伝子の同定を行います。
コマンドライン引数は3つ指定する必要があります。1つ目は、アノテーション情報（gtfファイル）を指定し、2つ目と3つ目には比較したいデータ（bamファイル）をそれぞれ指定します。Duplicateのデータ（n >= 2）がある場合は、各データ（bamファイル）を「,」カンマで区切り指定することができます（比較するデータをn = 1でしか取っていない場合、有意差を計算することができません。そのため、この場合Fold-change（例えば2以上）をもとに発現変動遺伝子を絞り込む必要があります）。
〈シェルスクリプト実行例〉
/path/to/RNA-seq_pipeline.sh SRRxxxxxx SRRyyyyyy SRRzzzzzz
⇒シェルスクリプトを呼び出し、コマンドライン引数として処理したいFastqファイルの名前を指定。
2013年6月10日月曜日
blastclustを利用した冗長配列のクラスタリング
BlastclustはNCBIが提供するBlastアルゴリズムを利用した塩基配列もしくはアミノ酸配列のクラスタリングに利用されるソフトです。
近年の論文でもBlastclustは活用されており、例えばENCODEプロジェクトにおいてncRNAのTranscriptomeデータをクラスタリングする際に用いられたソフトでもあります。
主な使用方法としては、クラスタリングしたい配列データをFASTAファイルの形で用意して、Blastclustを実行させるだけです。もちろん、目的に応じてパラメータをいじくる必要はあります。
〈Blastclustのダウンロード〉
（１）下記のURLからBlastパッケージをダウンロードする
ftp://ftp.ncbi.nih.gov/blast/executables/release/2.2.26/
⇒自分が使っているOS向けのパッケージを選択。
⇒近年、「Blast」の機能を向上させた「Blast+」がリリースされていますが、Blast+のパッケージには「Blastclust」が入っていません（もしかしたら、別名で入っているのかな？）。そのため、現時点で最も最新の「Blast 2.2.26」をここではダウンロード先としています。
〈Blastclustのインストール〉
Blastのソフトウェアはすでにビルドされているので、ファイルの解凍とパスを通すだけで使用できるようになります。
ただ、そのままファイルを解凍してしまうと、カレントディレクトリにBlastのbinディレクトリなどの中身が出力されてしまうので、あらかじめBlast用のディレクトリを作成し、その中でファイルを解凍することをおすすめします。
〈Blastclustの使用例〉
ここでは、遺伝子のisoformをクラスタリングしてみたいと思います。パラメータはENCODEプロジェクトで使用されたものをそのまま流用します。
〈参考文献〉
-blastclust manual page
http://www.csc.fi/english/research/sciences/bioscience/programs/blast/blastclust
-バイオインフォマティクス blastclust (legacy blast)
http://bi.biopapyrus.net/app/blast/blastclust.html
近年の論文でもBlastclustは活用されており、例えばENCODEプロジェクトにおいてncRNAのTranscriptomeデータをクラスタリングする際に用いられたソフトでもあります。
主な使用方法としては、クラスタリングしたい配列データをFASTAファイルの形で用意して、Blastclustを実行させるだけです。もちろん、目的に応じてパラメータをいじくる必要はあります。
〈Blastclustのダウンロード〉
（１）下記のURLからBlastパッケージをダウンロードする
ftp://ftp.ncbi.nih.gov/blast/executables/release/2.2.26/
⇒自分が使っているOS向けのパッケージを選択。
⇒近年、「Blast」の機能を向上させた「Blast+」がリリースされていますが、Blast+のパッケージには「Blastclust」が入っていません（もしかしたら、別名で入っているのかな？）。そのため、現時点で最も最新の「Blast 2.2.26」をここではダウンロード先としています。
〈Blastclustのインストール〉
Blastのソフトウェアはすでにビルドされているので、ファイルの解凍とパスを通すだけで使用できるようになります。
ただ、そのままファイルを解凍してしまうと、カレントディレクトリにBlastのbinディレクトリなどの中身が出力されてしまうので、あらかじめBlast用のディレクトリを作成し、その中でファイルを解凍することをおすすめします。
# BLASTディレクトリの作成 $ mkdir BLAST # ダウンロードしたblastパッケージをBLASTディレクトリに移動 $ mv ./blast-2.2.9-ia32-linux.tar.gz ./BLAST⇒「BLASTディレクトリ」を作成した親ディレクトリ内に「blast-2.2.9-ia32-linux.tar.gz」ファイルがあることを想定しています。
# BLASTディレクトリに移動 $ cd ./BLAST # ファイルの解凍 $ tar zxvf ./blast-2.2.9-ia32-linux.tar.gz
# ./profileのあるディレクトリに移動（Ubuntuのケース） $ cd $ pwd /home/imamachi # viエディタを用いた「./profile」の編集 $ vi .profile export PATH="$PATH:/path/to/BLAST/bin"を追加⇒geditなどのテキストエディタを用いて、編集してもよい。その場合、「.profile」は隠しファイルなので、Ubuntuのケースではホームフォルダからメニューバーの「表示」→「隠しファイルを表示する」にチェックを付けることで「.profile」を表示させる。
〈Blastclustの使用例〉
ここでは、遺伝子のisoformをクラスタリングしてみたいと思います。パラメータはENCODEプロジェクトで使用されたものをそのまま流用します。
$ blastclust -i /path/to/INPUT.fasta -o /path/to/OUTPUT.fasta \ -p F -L 0.4 -S 40 -i: 入力ファイルの指定（FASTA形式） -o: 出力ファイルの指定 -p: 塩基配列なら「F」、アミノ酸配列なら「P」 -L: アライメントする配列の長さのカバー率（0.0-1.0（0%-100%の意味）の数字を指定できる。例えば、40%なら「0.4」と記述する） -S:アライメントする配列の類似度（3以上の数字を指定した場合、identical residuesのパーセンテージを指定（3-100の範囲））⇒各パラメータの詳しい計算方法については、下記のサイトを参照。
〈参考文献〉
-blastclust manual page
http://www.csc.fi/english/research/sciences/bioscience/programs/blast/blastclust
-バイオインフォマティクス blastclust (legacy blast)
http://bi.biopapyrus.net/app/blast/blastclust.html
2013年5月26日日曜日
FASTX-toolkitをインストールする（Linux OS, Ubuntu）
FASTX-toolkitはFASTA・FASTQファイルの処理を行うことができるソフトウェアです。クオリティの高いリードを選抜したり、タグ配列の除去などに利用できます。
〈FASTX-toolkitのダウンロード〉
（１）下記のFASTX-toolkitのサイトへ移動。
http://hannonlab.cshl.edu/fastx_toolkit/index.html
（２）メニューの「Download & Installation」をクリック。
http://hannonlab.cshl.edu/fastx_toolkit/download.html
（３）「Download」の項にファイルが置いてあるので、
・fastx_toolkit-0.0.13.2.tar.bz2
・libgtextutils-0.6.1.tar.bz2
をそれぞれダウンロードする。
〈FASTX-toolkitのインストール〉
FASTX-toolkitには「libgtextutils-0.6」が必要となります。そこで、まずこのファイルをインストールしましょう。
⇒今回はもともとパスの通っている「/bin」ディレクトリにファイルをコピペしましたが、もちろん「.profile」ファイルにパスを通す先のディレクトリの場所を書き加えてもOKです。
FASTX-toolkitのコマンドのうち、
・fasta_clipping_histogram
・The fastx_barcode_splitter
・The fastq_quality_boxplot
の３つについては追加でPerlのモジュールやいくつかのソフトウェアをインストールする必要があります。個人的に、これら３つは今まで使ったことがないので、インストール方法については割愛します。
よくよくマニュアルをみてみると、「The fastx_barcode_splitter」を使ってFASTQファイルを分割できるようで、小さく分割したFASTQファイルを並列処理させる場合に使えそうですね。
〈参考文献〉
-FASTX-Toolkit_Command-line Usage
http://hannonlab.cshl.edu/fastx_toolkit/commandline.html
〈FASTX-toolkitのダウンロード〉
（１）下記のFASTX-toolkitのサイトへ移動。
http://hannonlab.cshl.edu/fastx_toolkit/index.html
（２）メニューの「Download & Installation」をクリック。
http://hannonlab.cshl.edu/fastx_toolkit/download.html
（３）「Download」の項にファイルが置いてあるので、
・fastx_toolkit-0.0.13.2.tar.bz2
・libgtextutils-0.6.1.tar.bz2
をそれぞれダウンロードする。
〈FASTX-toolkitのインストール〉
FASTX-toolkitには「libgtextutils-0.6」が必要となります。そこで、まずこのファイルをインストールしましょう。
#ファイルの解凍 $ bzip -d libgtextutils-0.6.1.tar.bz2 $ tar xvf libgtextutils-0.6.1.tar #ディレクトリの移動 $ cd ./libgtextutils-0.6.1 #コンパイル（インストール） $ ./configure $ make $ sudo make install
#ファイルの解凍 $ bzip -d libgtextutils-0.6.1 fastx_toolkit-0.0.13.2.tar.bz2 $ tar xvf libgtextutils-0.6.1 fastx_toolkit-0.0.13.2.tar #ディレクトリの移動 cd ./fastx_toolkit-0.0.13.2 #コンパイル $ ./configure $ make #パスを通す $sudo cp ./bin/* /bin⇒「*」は任意の文字列を表すメタ文字（ここでは、すべてのファイル名とマッチする）。
⇒今回はもともとパスの通っている「/bin」ディレクトリにファイルをコピペしましたが、もちろん「.profile」ファイルにパスを通す先のディレクトリの場所を書き加えてもOKです。
FASTX-toolkitのコマンドのうち、
・fasta_clipping_histogram
・The fastx_barcode_splitter
・The fastq_quality_boxplot
の３つについては追加でPerlのモジュールやいくつかのソフトウェアをインストールする必要があります。個人的に、これら３つは今まで使ったことがないので、インストール方法については割愛します。
よくよくマニュアルをみてみると、「The fastx_barcode_splitter」を使ってFASTQファイルを分割できるようで、小さく分割したFASTQファイルを並列処理させる場合に使えそうですね。
〈参考文献〉
-FASTX-Toolkit_Command-line Usage
http://hannonlab.cshl.edu/fastx_toolkit/commandline.html
ディスクの容量に困ったらファイルを圧縮しておこう
大規模データを処理していると、あっという間にディスクの容量を食います。そのため、どうにかしてディスクの容量を抑える必要があります。昔解析したデータを順次消していくというのが手っ取り早い方法かも知れませんが、過去に解析したデータを残しておきたいケースもあると思います。
そこで、以下の圧縮率が高い「bzip2」形式での圧縮方法を紹介。
ディレクトリごと圧縮・解凍する方法
j: bzip2形式の圧縮を同時に行う
v: 処理情報をprintする
f: アーカイブファイルを指定する
x: 書庫からファイルを取り出す
また、圧縮や解凍をしつつ他の作業をしたい場合、「nohup」コマンドを利用すると便利です。このコマンドを先頭につけることにより、端末を閉じても、Linux OS自体をシャットダウンしない限り、裏でコマンドでの処理をし続けてくれます。
nohupを使うとさらに便利
そこで、以下の圧縮率が高い「bzip2」形式での圧縮方法を紹介。
ディレクトリごと圧縮・解凍する方法
#圧縮 tar -cjvf dir.bz2.tar dirc: アーカイブファイルの作成
j: bzip2形式の圧縮を同時に行う
v: 処理情報をprintする
f: アーカイブファイルを指定する
#解凍 tar jxvf dir.bz2.tar-j: bzip2を通して処理を行う
x: 書庫からファイルを取り出す
また、圧縮や解凍をしつつ他の作業をしたい場合、「nohup」コマンドを利用すると便利です。このコマンドを先頭につけることにより、端末を閉じても、Linux OS自体をシャットダウンしない限り、裏でコマンドでの処理をし続けてくれます。
nohupを使うとさらに便利
nohup tar -cjvf dir.bz2.tar dir
2013年4月21日日曜日
Cytoscapeを利用したPathway解析をやってみた（Cytoscape×Reactome）
Cytoscapeとは、Gene ontology解析やPathway解析などのネットワーク構造を持つデータを可視化するソフトです。（主に、RNA-seqなどの発現量解析以降の下流の解析に利用可能なのではないかと思います。）
Cytoscapeにはさまざまなプラグインが存在し、プラグインを追加することにより様々な解析がCytoscape上で可能になります。今回は、「Reactome」のプラグインを導入して、（なんちゃって）Pathway解析をやってみたいと思います。
〈Cytoscapeのダウンロード〉
（１）Cytoscapeのサイトから上部のメニューの「Download」をクリック
（２）Your informationの「Name」「Organization」「Email address」の欄を埋めて、Conditionsの「I accept terms of use」にチェックを入れる。
（３）Download Latest Production Version: 2.8.3もしくはDownload New 3.x Series: 3.0.0の「Platform Specific installers」から自分のOSに合ったインストーラーを選択し、ファイルをダウンロードする。
〈Cytoscapeのインストール〉
（１）Windows版なら、「.exe」ファイルなのでNextボタンを適当にクリックしていくだけでインストール可能。
〈Reactomeのプラグインのダウンロード〉
（１）Cytoscape用プラグインを下記のURLからダウンロードする。
http://wiki.reactome.org/index.php/Reactome_FI_Cytoscape_Plugin
⇒caBigR3.jarをダウンロードするように文章中に書かれているのでこのファイルをダウンロードしてくる。
〈Reactomeのプラグインのインストール〉
（１）Windows版では、「C:/Program Files/Cytoscape_v2.8.3」にCytoscapeのフォルダが存在し、そのフォルダ内の「Plugins」フォルダの中に先ほどダウンロードした「caBigR3.jar」をコピペする（もしくは移動させる）。
（２）Cytoscapeを再起動させる。
〈Reactomeデータベースを利用したPathway解析〉
（１）Cytoscapeを立ち上げる。
（２）メニューバーから[Plugins]→[Reactome FIｓ]→[Gene Set/Mutation Analysis]を選択し、クリック。
（３）Loading Altered Genes中の「Browse」ボタンをクリックし、インポートしたい発現変動遺伝子のシンボルの書かれたファイルを選択する。
⇒発現変動遺伝子のリストの例は以下のサンプルファイルを参照。
1. Simple gene set: GWASFuzzyGenes.txt
2. Gene/sample number pair: GeneSampleNumber.txt
3. NCI MAF (mutation annotation file): GlioblastomaMutationTable.txt
（４）「OK」ボタンをクリックすると、ReactomeのPathwayデータベースを参照して、Node（節）とLine（線）で描かれた遺伝子の関係図を作ってくれる。
このReactome以外にも、さまざまなプラグインが存在し、論文も出ているので色々試してみようかと思います。中にはライセンス登録が必要なプラグインもあり、導入するのが面倒なものもあるようです。
〈参考文献〉
-Cytoscape: a software environment for integrated models of biomolecular interaction networks
http://www.ncbi.nlm.nih.gov/pubmed/14597658
-Cytoscape 2.8: new features for data integration and network visualization
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3031041/
-Reactome: a database of reactions, pathways and biological processes
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3013646/
アメリエフのブログ -オープンソースでパスウェイ解析（2013.04.17 Wednesday）
http://blog.amelieff.jp/
Cytoscapeにはさまざまなプラグインが存在し、プラグインを追加することにより様々な解析がCytoscape上で可能になります。今回は、「Reactome」のプラグインを導入して、（なんちゃって）Pathway解析をやってみたいと思います。
〈Cytoscapeのダウンロード〉
（１）Cytoscapeのサイトから上部のメニューの「Download」をクリック
（２）Your informationの「Name」「Organization」「Email address」の欄を埋めて、Conditionsの「I accept terms of use」にチェックを入れる。
（３）Download Latest Production Version: 2.8.3もしくはDownload New 3.x Series: 3.0.0の「Platform Specific installers」から自分のOSに合ったインストーラーを選択し、ファイルをダウンロードする。
〈Cytoscapeのインストール〉
（１）Windows版なら、「.exe」ファイルなのでNextボタンを適当にクリックしていくだけでインストール可能。
〈Reactomeのプラグインのダウンロード〉
（１）Cytoscape用プラグインを下記のURLからダウンロードする。
http://wiki.reactome.org/index.php/Reactome_FI_Cytoscape_Plugin
⇒caBigR3.jarをダウンロードするように文章中に書かれているのでこのファイルをダウンロードしてくる。
〈Reactomeのプラグインのインストール〉
（１）Windows版では、「C:/Program Files/Cytoscape_v2.8.3」にCytoscapeのフォルダが存在し、そのフォルダ内の「Plugins」フォルダの中に先ほどダウンロードした「caBigR3.jar」をコピペする（もしくは移動させる）。
（２）Cytoscapeを再起動させる。
〈Reactomeデータベースを利用したPathway解析〉
（１）Cytoscapeを立ち上げる。
（２）メニューバーから[Plugins]→[Reactome FIｓ]→[Gene Set/Mutation Analysis]を選択し、クリック。
（３）Loading Altered Genes中の「Browse」ボタンをクリックし、インポートしたい発現変動遺伝子のシンボルの書かれたファイルを選択する。
⇒発現変動遺伝子のリストの例は以下のサンプルファイルを参照。
1. Simple gene set: GWASFuzzyGenes.txt
2. Gene/sample number pair: GeneSampleNumber.txt
3. NCI MAF (mutation annotation file): GlioblastomaMutationTable.txt
（４）「OK」ボタンをクリックすると、ReactomeのPathwayデータベースを参照して、Node（節）とLine（線）で描かれた遺伝子の関係図を作ってくれる。
このReactome以外にも、さまざまなプラグインが存在し、論文も出ているので色々試してみようかと思います。中にはライセンス登録が必要なプラグインもあり、導入するのが面倒なものもあるようです。
〈参考文献〉
-Cytoscape: a software environment for integrated models of biomolecular interaction networks
http://www.ncbi.nlm.nih.gov/pubmed/14597658
-Cytoscape 2.8: new features for data integration and network visualization
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3031041/
-Reactome: a database of reactions, pathways and biological processes
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3013646/
アメリエフのブログ -オープンソースでパスウェイ解析（2013.04.17 Wednesday）
http://blog.amelieff.jp/
2013年4月20日土曜日
Motif Discovery（Weeder2weblogo）
Weederの「OUTPUTファイル」から、weblogoの「INPUTファイル」であるfastaファイルを生成するスクリプトを書いてみました。
Weederで予測されたモチーフ配列を視覚化するために、予測されたモチーフ配列に対してそれぞれのweblogo用のfastaファイルが生成させるスクリプトです。
このスクリプトに通した後、生成されたfastaファイルを一つ一つweblogoで視覚化。
・github -non-bioinformatician_NGS_analysis
https://github.com/Naoto-Imamachi/non-bioinformatician_NGS_analysis
Weederで予測されたモチーフ配列を視覚化するために、予測されたモチーフ配列に対してそれぞれのweblogo用のfastaファイルが生成させるスクリプトです。
このスクリプトに通した後、生成されたfastaファイルを一つ一つweblogoで視覚化。
#!/usr/local/bin/perl
####################################
#Weeder2weblogo#####################
#var1.0_Perl
#update: 2012.4.8
####################################
use strict;
use warnings;
my $file = $ARGV[0];
my $number = $ARGV[1];
#main####################################################
open (IN,"<./$file\.wee") || die;
my $flg_gene=0;
my $flg=0;
my $filehandle = 1;
my $count = 1;
my %gene_ID;
while(my $line = <IN>){
#Gene_ID
if($line =~ /^Your sequences/){
$flg_gene = 1;
next;
}elsif($flg_gene == 1){
$flg_gene = 2;
next;
}elsif($flg_gene == 2 && $line eq "\n"){
$flg_gene = 0;
next;
}elsif($flg_gene == 2){
my @data = split(/\s+/,$line);
$gene_ID{$data[1]}=$data[3];
next;
}
#Sequence_data
if($line =~ /^Best occurrences/){
if($filehandle == 0){
close(OUT);
$count++;
}else{
$filehandle = 0;
}
open (OUT,">./$file\#$count\.fasta") || die;
$flg = 1;
next;
}elsif($flg == 1){
$flg = 2;
next;
}elsif($line eq "\n" && $flg == 2){
$flg = 0;
next;
}elsif($flg == 2){
if($line =~ /^\s+/){
$line =~ s/^\s+//;
}
my @data = split(/\s+/,$line);
if($data[2] =~ /^\[/){
$data[2] =~ s/\[//;
$data[2] =~ s/\]//;
}
print OUT "$gene_ID{$data[0]}\n";
print OUT "$data[2]\n";
next;
}
}
close(IN);
close(OUT);
print "...process successfully completed\n";
〈使用方法〉
$ perl Weeder2weblogo INPUT_file ⇒INPUT_fileの拡張子は「.wee」にしておくこと。 （Weederからの出力ファイル「.wee」を利用しているので、拡張子を変更しないように。）〈ダウンロード先〉
・github -non-bioinformatician_NGS_analysis
https://github.com/Naoto-Imamachi/non-bioinformatician_NGS_analysis
2013年4月6日土曜日
BEDToolsをインストールする（Linux OS, Ubuntu）
「BEDTools」とは、SAM, GTF, BEDなどのファイル形式のデータを加工・処理するために開発された解析支援ツールです。
BEDツール同士を比較したり、BEDファイルからFASTAファイルへの変換を行ったり、BAMファイルからBEDファイルへの変換を行ったりできます。
いろいろな機能があるので、それらについての詳細はおいおい説明していければと思います。
今回は、その「BEDTools」のインストール方法についてのメモになります。
〈BEDToolsのダウンロード〉
（１）下記のbedtoolsのサイトへ移動。
https://code.google.com/p/bedtools/
（２）「Downloads」をクリック。
https://code.google.com/p/bedtools/downloads/list
（３）「BEDTools.v2.17.0.tar.gz」と書かれたファイルが置いてあるので、それをクリックする。
BEDTools.v2.17.0.tar.gz
（４）「File: BEDTools.v2.17.0.tar.gz 1.2MB」の部分をクリック。
https://bedtools.googlecode.com/files/BEDTools.v2.17.0.tar.gz
〈BEDToolsのインストール〉
BEDツール同士を比較したり、BEDファイルからFASTAファイルへの変換を行ったり、BAMファイルからBEDファイルへの変換を行ったりできます。
いろいろな機能があるので、それらについての詳細はおいおい説明していければと思います。
今回は、その「BEDTools」のインストール方法についてのメモになります。
〈BEDToolsのダウンロード〉
（１）下記のbedtoolsのサイトへ移動。
https://code.google.com/p/bedtools/
（２）「Downloads」をクリック。
https://code.google.com/p/bedtools/downloads/list
（３）「BEDTools.v2.17.0.tar.gz」と書かれたファイルが置いてあるので、それをクリックする。
BEDTools.v2.17.0.tar.gz
（４）「File: BEDTools.v2.17.0.tar.gz 1.2MB」の部分をクリック。
https://bedtools.googlecode.com/files/BEDTools.v2.17.0.tar.gz
〈BEDToolsのインストール〉
$ tar zxvf BEDTools.v2.17.0.tar.gz
$ cd bedtools-2.17.0/
$ make
#「.profile」ファイルを書き換える方法 $ vi ~/.profile ⇒最終行に、「export PATH="$PATH:/home/bowtie.XX.XX/…/任意のパス/bedtools-2.17.0"」を追加
$ bedtools --version bedtools v2.17.0 ⇒バージョンの確認。
LinuxとWindowsの改行コードの違いについて
LinuxとWindowsでは改行コードが違います。
Linuxでは、「ラインフィード（LF）」
Windowsでは、「キャリッジリターン＋ラインフィード（CR+LF）」
と呼ばれる改行コードがそれぞれ使われています。
テキストの中身をみてみると、
Linuxでは、「\n」
Windowsでは、「\r\n」
となっています。
このように、LinuxとWindowsでは改行コードが異なるために、Perlなどのスクリプト言語を使ったり、LinuxとWindows間でデータのやり取りしたりしていると、この改行コードの違いで問題を生じるケースがあります。
ちなみに、この「\」なんですが日本のキーボードだと二種類存在します(´・ω・`)
１つは、キーボード上右端の「￥キー」、
もう１つは、キーボード下右端の「＼キー」です。
Windows上では、「￥キー」も「＼キー」もPerlなどでは「\n」として認識されます。
しかし、Linux上ではどうやら「￥キー」＋「n」だと、改行コードとして認識してくれず、「￥n」という文字列として表示されてしまいます。
なので、Linux上でPerlなどのスクリプトを実行する際には、ちゃんと「＼キー」＋「n」と書く必要があります。
以下では、Linux上でのPerlスクリプトの実行例を示しています。
上記では、「￥」を全角で記述していますが、半角だとこのブログ上ではバックスラッシュとして表現されてしまうため、わざと全角にしてます。
Linuxでは、「ラインフィード（LF）」
Windowsでは、「キャリッジリターン＋ラインフィード（CR+LF）」
と呼ばれる改行コードがそれぞれ使われています。
テキストの中身をみてみると、
Linuxでは、「\n」
Windowsでは、「\r\n」
となっています。
このように、LinuxとWindowsでは改行コードが異なるために、Perlなどのスクリプト言語を使ったり、LinuxとWindows間でデータのやり取りしたりしていると、この改行コードの違いで問題を生じるケースがあります。
ちなみに、この「\」なんですが日本のキーボードだと二種類存在します(´・ω・`)
１つは、キーボード上右端の「￥キー」、
もう１つは、キーボード下右端の「＼キー」です。
Windows上では、「￥キー」も「＼キー」もPerlなどでは「\n」として認識されます。
しかし、Linux上ではどうやら「￥キー」＋「n」だと、改行コードとして認識してくれず、「￥n」という文字列として表示されてしまいます。
なので、Linux上でPerlなどのスクリプトを実行する際には、ちゃんと「＼キー」＋「n」と書く必要があります。
以下では、Linux上でのPerlスクリプトの実行例を示しています。
print "Hello￥n"; print "Have a nice day￥n"; # 実行 hello￥nhava a nice day￥n
print "Hello\n"; print "Have a nice day\n"; # 実行 hello hava a nice day※注意※
上記では、「￥」を全角で記述していますが、半角だとこのブログ上ではバックスラッシュとして表現されてしまうため、わざと全角にしてます。
パスの通し方2（Linux, Ubuntu）
パスの通し方の１つとしては「.profile」の最終行に「export PATH=$PATH:/任意のパス/:/任意のパス/:…」を書き加えるという方法があります。
この他にも、もっと簡単な方法として、パスの通っているディレクトリにコンパイル済みのプログラム（バイナリファイル）をコピペするという方法もあります。
以上。
この他にも、もっと簡単な方法として、パスの通っているディレクトリにコンパイル済みのプログラム（バイナリファイル）をコピペするという方法もあります。
$ sudo cp ./任意のプログラム名 /usr/local/bin [sudo] password for imamachi:⇒ここでパスワードを入力し、一時的に管理者権限を取得する。
以上。
Motif Discovery（BioProspector2weblogo）
BioProspectorの「OUTPUTファイル」から、weblogoの「INPUTファイル」であるfastaファイルを生成するスクリプトを書いてみました。
BioProspectorでは、複数のモチーフ配列が予想されます。このスクリプトに通すことで予測されたモチーフ配列に対してそれぞれのweblogo用のfastaファイルが生成されるようになっています。（５つのモチーフ配列が予想されていれば、このスクリプトで５つのfastaファイルが生成されます。）
あとは、生成されたfastaファイルを一つ一つweblogoで視覚化するだけ。（あと、一連のステップもバッチ処理で一気に生成出来ればよいのですが…。）
〈ダウンロード先〉
・github -non-bioinformatician_NGS_analysis
https://github.com/Naoto-Imamachi/non-bioinformatician_NGS_analysis
BioProspectorでは、複数のモチーフ配列が予想されます。このスクリプトに通すことで予測されたモチーフ配列に対してそれぞれのweblogo用のfastaファイルが生成されるようになっています。（５つのモチーフ配列が予想されていれば、このスクリプトで５つのfastaファイルが生成されます。）
あとは、生成されたfastaファイルを一つ一つweblogoで視覚化するだけ。（あと、一連のステップもバッチ処理で一気に生成出来ればよいのですが…。）
#!/usr/local/bin/perl
####################################
#BioProspector2weblogo##############
#var1.0_Perl
#update: 2012.4.5
####################################
use strict;
use warnings;
my $file = $ARGV[0];
my $number = $ARGV[1];
#main####################################################
open (IN,"<./$file\.result") || die;
my $flg=0;
my $filehandle = 1;
my $count = 1;
while(my $line = <in>){
if($line eq "\n"){
next;
}
chomp $line;
if($line =~ /^Try/
|| $line =~ /^The highest/
|| $line =~ /^Width/
|| $line =~ /^Blk1/
|| $line =~ /^[1-6]/){
next;
}
if($line =~ /^Motif/){
if($filehandle == 0){
close(OUT);
$count++;
}else{
$filehandle = 0;
}
open (OUT,">./$file\#$count\.fasta") || die;
}else{
if($line =~ /^>/){
print OUT "$line\n";
$flg=1;
}elsif($flg == 1){
print OUT "$line\n";
$flg=0;
}
}
}
close(IN);
close(OUT);
print "...process successfully completed\n";
〈使用方法〉$ perl Bioprospector2weblogo INPUT_file⇒INPUT_fileの拡張子は「.result」にしておくこと。（もしくはスクリプトを適宜書き換え。）
〈ダウンロード先〉
・github -non-bioinformatician_NGS_analysis
https://github.com/Naoto-Imamachi/non-bioinformatician_NGS_analysis
2013年3月31日日曜日
Sequence Logoの作成（WebLogoを使ってみた）
シークエンスロゴ（Sequence Logo）とは、特定の部位における塩基もしくはアミノ酸の保存性（もしくは出現頻度）を視覚化した図のことです。
シークエンスロゴによるモチーフ配列（コンセンサス配列）の表示に使われたりします。
〈WebLogoに必要なパッケージのインストール〉
Weblogoをインストールする前に、「ghostscript」「pdf2svg」「numpy」とpythonの「setuptools」をあらかじめインストールする必要があります。以下、「apt-get」commandを利用してパッケージをインストールしていきたいと思います。
（１）ghostscriptのインストール
$ apt-file search /usr/bin/ghostscript
ghostscript: /usr/bin/ghostscript
⇒パッケージ名を検索。
$ sudo apt-get install ghostscript
[sudo] password for imamachi:
⇒パッケージのインストール。管理者権限により実行。
⇒Ubuntuではデフォルトでインストールされている場合があります。
（２）pdf2svgのインストール
$ apt-file search pdf2svg
docvert-libreoffice: /usr/share/docvert/core/config/unix-specific/convert-using-pdf2svg.sh
pdf2svg: /usr/bin/pdf2svg
pdf2svg: /usr/share/doc/pdf2svg/changelog.Debian.gz
pdf2svg: /usr/share/doc/pdf2svg/copyright
pdf2svg: /usr/share/man/man1/pdf2svg.1.gz
ruby-poppler: /usr/share/doc/ruby-poppler/examples/pdf2svg.rb
$ sudo apt-get install pdf2svg
[sudo] password for imamachi:
（３）numpyのインストール
$ apt-file search usr/include/numpy
python-numpy: /usr/include/numpy
$ sudo apt-get install python-numpy
[sudo] password for imamachi:
（４）pythonのsetuptoolsのインストール
$ apt-file search setuptools
python-setuptools: /usr/lib/python2.7/dist-packages/setuptools.egg-info
python-setuptools: /usr/lib/python2.7/dist-packages/setuptools.pth
python-setuptools: /usr/lib/python2.7/dist-packages/setuptools/__init__.py
python-setuptools: /usr/lib/python2.7/dist-packages/setuptools/archive_util.py
python-setuptools: /usr/lib/python2.7/dist-packages/setuptools/cli-32.exe
…
$ sudo apt-get install python-setuptools
[sudo] password for imamachi:
〈WebLogoのインストール〉
（１）WebLogoをpythonのsetuptoolsを利用し、インストールする。
$ sudo easy_install weblogo
[sudo] password for imamachi:
Searching for weblogo
Reading http://pypi.python.org/simple/weblogo/
Reading http://code.google.com/p/weblogo/
Best match: weblogo 3.3
Downloading http://weblogo.googlecode.com/files/weblogo-3.3.tar.gz
Processing weblogo-3.3.tar.gz
Running weblogo-3.3/setup.py -q bdist_egg --dist-dir /tmp/easy_install-uW4jWu/weblogo-3.3/egg-dist-tmp-qHfZq3
warning: no files found matching '*.html' under directory 'weblogolib/htdocs/'
warning: no files found matching '*.cgi' under directory 'weblogolib/htdocs/'
warning: no files found matching '*.*' under directory 'weblogolib/htdocs/'
warning: no files found matching '*.*' under directory 'test_weblogo/data/'
warning: no files found matching '*.*' under directory 'test_corebio/'
warning: no files found matching '*.*' under directory 'test_corebio/data/'
warning: no files found matching '*.html' under directory 'apidocs/'
warning: no files found matching '*' under directory 'corebio/data/'
warning: no files found matching '*.*' under directory 'corebio/data/'
zip_safe flag not set; analyzing archive contents...
corebio.data: module references __file__
weblogolib._cli: module references __file__
weblogolib.__init__: module references __file__
weblogolib._cgi: module references __file__
Adding weblogo 3.3 to easy-install.pth file
Installing weblogo script to /usr/local/bin
Installing transformseq script to /usr/local/bin
Installed /usr/local/lib/python2.7/dist-packages/weblogo-3.3-py2.7.egg
Processing dependencies for weblogo
Finished processing dependencies for weblogo
⇒途中、警告文が出ているが、インストール自体はうまく行っているよう。
（２）WebLogoがインストールされているか確認する。
$ weblogo --version
WebLogo 3.3 (2012-07-02)
〈Weblogoを実行する〉
$ weblogo -f INPUT_file.fasta -o OUTPUT_file.png -F png
⇒実行例。他のオプションについてはWebLogoのサイトを確認。
・オプション
-f: INPUTファイルの指定（fasta format）
-o: OUTPUTファイルの指定（とりあえず拡張子もつけておく。）
-F: OUTPUTファイルのファイルフォーマットの指定。
〈参考文献〉
・Sequence logo - Wikipedia
http://en.wikipedia.org/wiki/Sequence_logo
・WebLogo 3 : User's Manual
http://weblogo.threeplusone.com/manual.html
シークエンスロゴによるモチーフ配列（コンセンサス配列）の表示に使われたりします。
〈WebLogoに必要なパッケージのインストール〉
Weblogoをインストールする前に、「ghostscript」「pdf2svg」「numpy」とpythonの「setuptools」をあらかじめインストールする必要があります。以下、「apt-get」commandを利用してパッケージをインストールしていきたいと思います。
（１）ghostscriptのインストール
$ apt-file search /usr/bin/ghostscript
ghostscript: /usr/bin/ghostscript
⇒パッケージ名を検索。
$ sudo apt-get install ghostscript
[sudo] password for imamachi:
⇒パッケージのインストール。管理者権限により実行。
⇒Ubuntuではデフォルトでインストールされている場合があります。
（２）pdf2svgのインストール
$ apt-file search pdf2svg
docvert-libreoffice: /usr/share/docvert/core/config/unix-specific/convert-using-pdf2svg.sh
pdf2svg: /usr/bin/pdf2svg
pdf2svg: /usr/share/doc/pdf2svg/changelog.Debian.gz
pdf2svg: /usr/share/doc/pdf2svg/copyright
pdf2svg: /usr/share/man/man1/pdf2svg.1.gz
ruby-poppler: /usr/share/doc/ruby-poppler/examples/pdf2svg.rb
$ sudo apt-get install pdf2svg
[sudo] password for imamachi:
（３）numpyのインストール
$ apt-file search usr/include/numpy
python-numpy: /usr/include/numpy
$ sudo apt-get install python-numpy
[sudo] password for imamachi:
（４）pythonのsetuptoolsのインストール
$ apt-file search setuptools
python-setuptools: /usr/lib/python2.7/dist-packages/setuptools.egg-info
python-setuptools: /usr/lib/python2.7/dist-packages/setuptools.pth
python-setuptools: /usr/lib/python2.7/dist-packages/setuptools/__init__.py
python-setuptools: /usr/lib/python2.7/dist-packages/setuptools/archive_util.py
python-setuptools: /usr/lib/python2.7/dist-packages/setuptools/cli-32.exe
…
$ sudo apt-get install python-setuptools
[sudo] password for imamachi:
〈WebLogoのインストール〉
（１）WebLogoをpythonのsetuptoolsを利用し、インストールする。
$ sudo easy_install weblogo
[sudo] password for imamachi:
Searching for weblogo
Reading http://pypi.python.org/simple/weblogo/
Reading http://code.google.com/p/weblogo/
Best match: weblogo 3.3
Downloading http://weblogo.googlecode.com/files/weblogo-3.3.tar.gz
Processing weblogo-3.3.tar.gz
Running weblogo-3.3/setup.py -q bdist_egg --dist-dir /tmp/easy_install-uW4jWu/weblogo-3.3/egg-dist-tmp-qHfZq3
warning: no files found matching '*.html' under directory 'weblogolib/htdocs/'
warning: no files found matching '*.cgi' under directory 'weblogolib/htdocs/'
warning: no files found matching '*.*' under directory 'weblogolib/htdocs/'
warning: no files found matching '*.*' under directory 'test_weblogo/data/'
warning: no files found matching '*.*' under directory 'test_corebio/'
warning: no files found matching '*.*' under directory 'test_corebio/data/'
warning: no files found matching '*.html' under directory 'apidocs/'
warning: no files found matching '*' under directory 'corebio/data/'
warning: no files found matching '*.*' under directory 'corebio/data/'
zip_safe flag not set; analyzing archive contents...
corebio.data: module references __file__
weblogolib._cli: module references __file__
weblogolib.__init__: module references __file__
weblogolib._cgi: module references __file__
Adding weblogo 3.3 to easy-install.pth file
Installing weblogo script to /usr/local/bin
Installing transformseq script to /usr/local/bin
Installed /usr/local/lib/python2.7/dist-packages/weblogo-3.3-py2.7.egg
Processing dependencies for weblogo
Finished processing dependencies for weblogo
⇒途中、警告文が出ているが、インストール自体はうまく行っているよう。
（２）WebLogoがインストールされているか確認する。
$ weblogo --version
WebLogo 3.3 (2012-07-02)
〈Weblogoを実行する〉
$ weblogo -f INPUT_file.fasta -o OUTPUT_file.png -F png
⇒実行例。他のオプションについてはWebLogoのサイトを確認。
・オプション
-f: INPUTファイルの指定（fasta format）
-o: OUTPUTファイルの指定（とりあえず拡張子もつけておく。）
-F: OUTPUTファイルのファイルフォーマットの指定。
〈参考文献〉
・Sequence logo - Wikipedia
http://en.wikipedia.org/wiki/Sequence_logo
・WebLogo 3 : User's Manual
http://weblogo.threeplusone.com/manual.html
Command Line Interfaceからパッケージをインストール（apt-getの活用）
今回は、Synapticパッケージマネジャーを利用しないやり方の紹介になります。
何かのパッケージをLinuxにインストールするのに、いちいちウェブサイトに出向いて、目的のパッケージをダウンロードしてきて、インストールするのは億劫ですよね。
そこで、「apt-get」コマンドの出番です( ･`ω･´)
このコマンドで、目的のパッケージを自動的にダウンロード・インストールまで行なってくれます。また、足りないパッケージを追加でインストールしてくれるので、さらに便利。
〈apt-fileで目的のファイルを含むパッケージの候補を検索〉
（１）まず、インストールしたいファイル・パッケージを検索します。
$apt-file search 任意のファイル名
libexpat1-dev: /usr/include/expat.h #libexpat1-devがパッケージ名になります。
⇒様々の候補が出力されるので、任意のものを選択。
〈apt-getで目的のファイルを含むパッケージをインストール〉
（１）管理者権限のもとで、インストールを実行します。
$ sudo apt-get install 任意のパッケージ名
あとは、インストールが終わるのをぼーっと待っています。
使用するディスク容量が大きい場合は以下のような文章が表示されます。
この操作後に追加で 17.6 MB のディスク容量が消費されます。
続行しますか [Y/n]? y
パッケージをインストールしていいか聞いてくるので、「y」と入力して、「Enterキー」を叩きましょう。問題が発生しなければ、インストールが実行されます。
何かのパッケージをLinuxにインストールするのに、いちいちウェブサイトに出向いて、目的のパッケージをダウンロードしてきて、インストールするのは億劫ですよね。
そこで、「apt-get」コマンドの出番です( ･`ω･´)
このコマンドで、目的のパッケージを自動的にダウンロード・インストールまで行なってくれます。また、足りないパッケージを追加でインストールしてくれるので、さらに便利。
〈apt-fileで目的のファイルを含むパッケージの候補を検索〉
（１）まず、インストールしたいファイル・パッケージを検索します。
$apt-file search 任意のファイル名
libexpat1-dev: /usr/include/expat.h #libexpat1-devがパッケージ名になります。
⇒様々の候補が出力されるので、任意のものを選択。
〈apt-getで目的のファイルを含むパッケージをインストール〉
（１）管理者権限のもとで、インストールを実行します。
$ sudo apt-get install 任意のパッケージ名
あとは、インストールが終わるのをぼーっと待っています。
使用するディスク容量が大きい場合は以下のような文章が表示されます。
この操作後に追加で 17.6 MB のディスク容量が消費されます。
続行しますか [Y/n]? y
パッケージをインストールしていいか聞いてくるので、「y」と入力して、「Enterキー」を叩きましょう。問題が発生しなければ、インストールが実行されます。
端末内でのコピー＆ペーストを可能にする（Linux, Ubuntu, VMware）
Ubuntuの端末を使っていて困るのが、「CTRL+C」でコピー、「CTRL+V」で貼り付けができないことです。
そこで、Ubuntuの端末のキーボードショートカットキーの変更を行うことで、上記のショートカットキーを使えるように設定を変更したいと思います。
（１）まず、端末を開く。
（２）メニューバーから、[編集]→[キーボードショートカット]をクリック。
（３）ウインドウが開くので、操作欄の「編集」にある「コピー」と「貼り付け」のショートカットキーを任意のキーに変更する。
⇒項目を選択し、今のショートカットキーをクリックする。すると、「新しいアクセラレータ…」と表示されるので、その状態で任意のキーを入力することで、ショートカットキーが変更される。
そこで、Ubuntuの端末のキーボードショートカットキーの変更を行うことで、上記のショートカットキーを使えるように設定を変更したいと思います。
（１）まず、端末を開く。
（２）メニューバーから、[編集]→[キーボードショートカット]をクリック。
（３）ウインドウが開くので、操作欄の「編集」にある「コピー」と「貼り付け」のショートカットキーを任意のキーに変更する。
⇒項目を選択し、今のショートカットキーをクリックする。すると、「新しいアクセラレータ…」と表示されるので、その状態で任意のキーを入力することで、ショートカットキーが変更される。
fastqファイルのクオリティチェック（fastq_illumina_filterを使ってみた）
illuminaのCASAVA-1.8 pipeline（Version 1.8）で生成されたfastqファイルは、クオリティの低いリードには「Y」、クオリティの良いリードには「N」が与えられています。（クオリティの良し悪しはCASAVA filteringにより落ちたリードかどうかで判断。フィルタリングの具体的な方法はよくわかっていない。）
fastqファイルは１つの配列の情報が４行で記述されており、CASAVA filterから落ちたリードかどうかはその１行目に書かれています。（fastqファイルのデータの読み方についてはwikiを参照。）
NCBIの「Gene Expression Omnibus（GEO）」か、DDBJの「Sequence Read Archive（SRA）」からダウンロードしたsraファイルを「SRA Toolskit」でfastqファイルに変換すると、illuminaのCASAVA pipelineから生成されたfastqファイルであっても、CASAVA filterの情報（NもしくはY）が抜け落ちています。（fastq-dumpコマンドを実行すると、自動的にクオリティの低いリードはフィルタリングされているよう。そのため、あまり気にする必要がないのかも。）
そのため、GEOやSRAのような「Public database」からダウンロードしてきたsraファイル（fastqファイル）については、以下の作業を行う必要はないと思います。
ただ、illuminaのCASAVA-1.8 pipelineを使って自前でfastqファイルを生成した場合では、以下のCASAVA filteringでクオリティが低いと判断されたリード（「Y」と記載）を排除したほうが良いと思います。
〈fastq_illumina_filter-0.1のインストール〉
（１）wgetコマンドで指定したサイトからファイルをダウンロードする。
$ wget http://cancan.cshl.edu/labmembers/gordon/fastq_illumina_filter/fastq_illumina_filter-0.1.tar.gz
⇒もしくは以下のサイトから、「Download」の「Source code」にあるファイルをダウンロード。
http://cancan.cshl.edu/labmembers/gordon/fastq_illumina_filter/
（２）ファイルを解凍する。
$ tar -xzf fastq_illumina_filter-0.1.tar.gz
（３）解凍したファイルのディレクトリに移動する。
$ cd fastq_illumina_filter-0.1
（４）コンパイルを実行する。
$ make
g++ -O2 -g -Wall -Wextra -Werror fastq_illumina_filter.cpp -o fastq_illumina_filter
（５）パスの通っているディレクトリへ、プログラム（fastq_illumina_filter）を移動する。
$ sudo cp fastq_illumina_filter /usr/local/bin
[sudo] password for imamachi:
⇒一時的に管理者権限を利用する。
（６）インストールされているか、実際にコマンドを実行してみる。
$ fastq_illumina_filter
fastq_illumina_filter (version 0.1) - Filters a FASTQ file generated by CASAVA 1.8
Copyright (C) 2011 - A. Gordon <gordon@cshl.edu> - Released under AGPLv3 or later.
Usage: fastq_illumina_filter [--keep Y/N] [-NYhv] [-o OUTPUT] [INPUT]
[INPUT] = Input file. Reads from STDIN if no file is specified.
[-o OUTPUT] = Output file. Default is STDOUT if no file is specified.
--keep N = Keep reads that were NOT filtered.
(Reads that have 'N' in the read-ID line.)
--keep Y = Keep reads that were filtered.
(Reads that have 'Y' in the read-ID line.)
-N = same as '--keep N'
-Y = save as '--keep Y'
-v = Report read counts to STDERR
Use twice to show progress while procesing the file.
When combined with '-o FILE', report goes to STDOUT.
-h = This helpful help screen
!! NOTE ABOUT Y/N !!:
Reads that were filtered (have 'Y' in the read-ID)
are the LOW QUALITY reads. You most likely DO NOT want them.
(In previous CASAVA version, those are reads that have not passed filtering)
Reads that were NOT filtered are the better-quality reads.
So using '-N' or '--keep N' is probably the option you want to use.
FASTQ files from CASAVA-1.8 Should have the following READ-ID format:
@<instrument>:<run number>:<flowcell ID>:<lane>:<tile>:<x-pos>:<y-pos> <read>:<is filtered>:<control number>:<index sequence>
Example 1:
$ gunzip -c NA10831_ATCACG_L002_R1_001.fastq.gz | head -n 4
@EAS139:136:FC706VJ:2:5:1000:12850 1:N:18:ATCACG
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
BBBBCCCC?<A?BC?7@@???????DBBA@@@@A@@
$ gunzip -c NA10831_ATCACG_L002_R1_001.fastq.gz | fastq_illumina_filter -vvN | gzip > good_reads.fq.gz
Processed 4,000,000 reads
fastq_illumina_filter (--keep N) statistics:
Input: 4,000,000 reads
Output: 3,845,435 reads (96%)
Example 2 (The input and output files are not compressed):
$ fastq_illumina_filter --keep N -v -v -o good_reads.fq NA10831_ATCACG_L002_R1_001.fastq
Processed 4,000,000 reads
fastq_illumina_filter (--keep N) statistics:
Input: 4,000,000 reads
Output: 3,845,435 reads (96%)
〈使用例〉
$ fastq_illumina_filter -N -v -o OUTPUT_file.fastq INPUT_file.fastq
・オプション
-N: CASAVA filterで落ちなかったリードをOUTPUT_file.fastqに出力。
-v: リード数の情報などを画面に出力。（STDERR）
など
〈参考文献〉
・Fastq - Wikipedia
http://ja.wikipedia.org/wiki/Fastq
fastqファイルは１つの配列の情報が４行で記述されており、CASAVA filterから落ちたリードかどうかはその１行目に書かれています。（fastqファイルのデータの読み方についてはwikiを参照。）
NCBIの「Gene Expression Omnibus（GEO）」か、DDBJの「Sequence Read Archive（SRA）」からダウンロードしたsraファイルを「SRA Toolskit」でfastqファイルに変換すると、illuminaのCASAVA pipelineから生成されたfastqファイルであっても、CASAVA filterの情報（NもしくはY）が抜け落ちています。（fastq-dumpコマンドを実行すると、自動的にクオリティの低いリードはフィルタリングされているよう。そのため、あまり気にする必要がないのかも。）
そのため、GEOやSRAのような「Public database」からダウンロードしてきたsraファイル（fastqファイル）については、以下の作業を行う必要はないと思います。
ただ、illuminaのCASAVA-1.8 pipelineを使って自前でfastqファイルを生成した場合では、以下のCASAVA filteringでクオリティが低いと判断されたリード（「Y」と記載）を排除したほうが良いと思います。
〈fastq_illumina_filter-0.1のインストール〉
（１）wgetコマンドで指定したサイトからファイルをダウンロードする。
$ wget http://cancan.cshl.edu/labmembers/gordon/fastq_illumina_filter/fastq_illumina_filter-0.1.tar.gz
⇒もしくは以下のサイトから、「Download」の「Source code」にあるファイルをダウンロード。
http://cancan.cshl.edu/labmembers/gordon/fastq_illumina_filter/
（２）ファイルを解凍する。
$ tar -xzf fastq_illumina_filter-0.1.tar.gz
（３）解凍したファイルのディレクトリに移動する。
$ cd fastq_illumina_filter-0.1
（４）コンパイルを実行する。
$ make
g++ -O2 -g -Wall -Wextra -Werror fastq_illumina_filter.cpp -o fastq_illumina_filter
（５）パスの通っているディレクトリへ、プログラム（fastq_illumina_filter）を移動する。
$ sudo cp fastq_illumina_filter /usr/local/bin
[sudo] password for imamachi:
⇒一時的に管理者権限を利用する。
（６）インストールされているか、実際にコマンドを実行してみる。
$ fastq_illumina_filter
fastq_illumina_filter (version 0.1) - Filters a FASTQ file generated by CASAVA 1.8
Copyright (C) 2011 - A. Gordon <gordon@cshl.edu> - Released under AGPLv3 or later.
Usage: fastq_illumina_filter [--keep Y/N] [-NYhv] [-o OUTPUT] [INPUT]
[INPUT] = Input file. Reads from STDIN if no file is specified.
[-o OUTPUT] = Output file. Default is STDOUT if no file is specified.
--keep N = Keep reads that were NOT filtered.
(Reads that have 'N' in the read-ID line.)
--keep Y = Keep reads that were filtered.
(Reads that have 'Y' in the read-ID line.)
-N = same as '--keep N'
-Y = save as '--keep Y'
-v = Report read counts to STDERR
Use twice to show progress while procesing the file.
When combined with '-o FILE', report goes to STDOUT.
-h = This helpful help screen
!! NOTE ABOUT Y/N !!:
Reads that were filtered (have 'Y' in the read-ID)
are the LOW QUALITY reads. You most likely DO NOT want them.
(In previous CASAVA version, those are reads that have not passed filtering)
Reads that were NOT filtered are the better-quality reads.
So using '-N' or '--keep N' is probably the option you want to use.
FASTQ files from CASAVA-1.8 Should have the following READ-ID format:
@<instrument>:<run number>:<flowcell ID>:<lane>:<tile>:<x-pos>:<y-pos> <read>:<is filtered>:<control number>:<index sequence>
Example 1:
$ gunzip -c NA10831_ATCACG_L002_R1_001.fastq.gz | head -n 4
@EAS139:136:FC706VJ:2:5:1000:12850 1:N:18:ATCACG
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
+
BBBBCCCC?<A?BC?7@@???????DBBA@@@@A@@
$ gunzip -c NA10831_ATCACG_L002_R1_001.fastq.gz | fastq_illumina_filter -vvN | gzip > good_reads.fq.gz
Processed 4,000,000 reads
fastq_illumina_filter (--keep N) statistics:
Input: 4,000,000 reads
Output: 3,845,435 reads (96%)
Example 2 (The input and output files are not compressed):
$ fastq_illumina_filter --keep N -v -v -o good_reads.fq NA10831_ATCACG_L002_R1_001.fastq
Processed 4,000,000 reads
fastq_illumina_filter (--keep N) statistics:
Input: 4,000,000 reads
Output: 3,845,435 reads (96%)
〈使用例〉
$ fastq_illumina_filter -N -v -o OUTPUT_file.fastq INPUT_file.fastq
・オプション
-N: CASAVA filterで落ちなかったリードをOUTPUT_file.fastqに出力。
-v: リード数の情報などを画面に出力。（STDERR）
など
〈参考文献〉
・Fastq - Wikipedia
http://ja.wikipedia.org/wiki/Fastq
2013年3月30日土曜日
一秒ごとにメモリの使用量を確認する方法（Linux）
解析をしていて気がかりなのが、メモリーの使用量です。
特に、マッピングの際に使用されるメモリーをモニタリングして、メモリー不足になっていないか確認したいと思うことがあります。（しょぼいスペックしか持たないパソコンなら尚更です…。）
メモリの使用量を見るには「free」コマンド、プログラムを一定間隔で実行し続ける「watch」コマンドを利用します。コマンドは以下のとおり。
（１）一秒ごとにメモリの使用量をMB単位で表示させる。
$ watch -n 1 free -m
-n: 秒単位でコマンドを実行する間隔を指定する。
-m: MB単位でメモリーの使用量を確認する。
以上。
特に、マッピングの際に使用されるメモリーをモニタリングして、メモリー不足になっていないか確認したいと思うことがあります。（しょぼいスペックしか持たないパソコンなら尚更です…。）
メモリの使用量を見るには「free」コマンド、プログラムを一定間隔で実行し続ける「watch」コマンドを利用します。コマンドは以下のとおり。
（１）一秒ごとにメモリの使用量をMB単位で表示させる。
$ watch -n 1 free -m
-n: 秒単位でコマンドを実行する間隔を指定する。
-m: MB単位でメモリーの使用量を確認する。
以上。
2013年3月29日金曜日
Motif Discovery （Weederを使ってみた）
前回に引き続き、コンセンサス配列（モチーフ配列）の検索ツールについての説明です。
今回は、「Weeder」と呼ばれるソフトを使ってみたいと思います。このソフトもよく利用されていると思います。
〈インストール方法〉
（１）下記のサイトからWeederをダウンロードする。
http://159.149.160.51/modtools/
⇒左側の「Downloads」から「Weeder1.4.2」をクリックする。ダウンロードのページに移動するので、Unix/Linux/Cygwin Source Code (NEW: version 1.4.2)をクリックし、ファイルをダウンロードする。
（２）ファイルを解凍する。
$ tar zxvf weeder1.4.2.tar.gz
（３）Weederのディレクトリに移動する。
$ cd Weeder1.4.2/
（４）コンパイルを実行。
$ ./compileall
gcc src/weederTFBS.c -o weederTFBS.out -O2 -lm
src/weederTFBS.c: 関数 ‘main’ 内:
src/weederTFBS.c:1050:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 3 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:1050:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 4 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:1052:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 3 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:1052:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 4 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:1054:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 3 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:1054:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 4 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:458:10: 警告: warn_unused_result 属性付きで宣言されている ‘fscanf’ の戻り値を無視しています [-Wunused-result]
gcc src/weederlauncher.c -o weederlauncher.out -lm
gcc src/adviser.c -o adviser.out -O2 -lm
src/adviser.c: 関数 ‘main’ 内:
src/adviser.c:289:10: 警告: warn_unused_result 属性付きで宣言されている ‘fscanf’ の戻り値を無視しています [-Wunused-result]
gcc src/locator.c -o locator.out -O2 -lm
src/locator.c: 関数 ‘main’ 内:
src/locator.c:329:10: 警告: warn_unused_result 属性付きで宣言されている ‘fscanf’ の戻り値を無視しています [-Wunused-result]
⇒よくわからない警告文が出ていましたが、インストールはうまくいったようです。
〈Frequency fileの作成〉
Weederの場合、バックグラウンドとなるデータをあらかじめ作っておく必要があります。（転写因子の結合サイトを探す場合は、すでに用意されているファイルを使えばよいみたいです。）
（１）下記のサイトへもう一度行き、左側の「Additional tools」のFrequency maker (for adding new species to Weeder)をクリック。
http://159.149.160.51/modtools/
（２）「How to build a new frequency file for Weeder」というページに移動するので、上から２つ目の「Download this program (works under UNIX/LINUX/MacOSX/Cygwin)」を右クリック、リンク先のファイルを保存してください。（名前は「wfrequency_maker.cpp」のままでダウンロードする。）
（３）ファイルをダウンロードしたディレクトリに移動し、コンパイルを実行。
$ g++ wfrequency_maker.cpp -o wfrequency_maker
（４）バックグラウンドとして使用する配列のFASTAファイルのパス、およびバックグラウンドデータを呼び出す際に用いる任意の２文字を入力して、Frequency fileを作成する。
$ ./wfrequency_maker ./任意のパス/INPUT_file.fasta XX（任意の２文字）
Reading sequences file (avoiding duplicates)... done - Read sequence(s): 9360 Good Sequences: 9325 Duplicated sequence(s): 34
Building 6mers list... done
Finding exact occurrences... done
Counting 6mers occurrences...4096/4096 done
Writing freq file: UP.6.freq done
Building 8mers list... done
Finding exact occurrences... done
Counting 8mers occurrences...65527/65527 done
Writing freq file: UP.8.freq done
XX.6.freq
XX.8.freq
⇒上記の２種類のファイルができるので、これらのファイルをWeeder1.4.2ディレクトリ中の「FreqFiles」へ移動させる。
〈モチーフ検索〉
$ ./weederlauncher.out ./任意のパス/INPUT_file.fasta XX（バックグラウンドデータ） analysis_mode
⇒上記はあくまでコマンドの入力例になります。詳しくは、「weedermanual.pdf」を呼んでみてください。
・パラメータ
analysis_modeはいくつか用意されています。
small: 6mer、8merのモチーフを検索する
median: 6mer、8mer、10merのモチーフを検索する
large: 6mer、8mer、10mer、12merのモチーフを検索する
などなど。
・オプションパラメータ
M: モチーフ配列が同一の配列に対して「２回以上」出現するモチーフ配列を検索する。
T[1-n]: モチーフ検索の回数を指定する。（1-n回: default 10）
などなど。
〈参考文献〉
Weeder Web: discovery of transcription factor binding sites in a set of sequences from co-regulated genes
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC441603/
今回は、「Weeder」と呼ばれるソフトを使ってみたいと思います。このソフトもよく利用されていると思います。
〈インストール方法〉
（１）下記のサイトからWeederをダウンロードする。
http://159.149.160.51/modtools/
⇒左側の「Downloads」から「Weeder1.4.2」をクリックする。ダウンロードのページに移動するので、Unix/Linux/Cygwin Source Code (NEW: version 1.4.2)をクリックし、ファイルをダウンロードする。
（２）ファイルを解凍する。
$ tar zxvf weeder1.4.2.tar.gz
（３）Weederのディレクトリに移動する。
$ cd Weeder1.4.2/
（４）コンパイルを実行。
$ ./compileall
gcc src/weederTFBS.c -o weederTFBS.out -O2 -lm
src/weederTFBS.c: 関数 ‘main’ 内:
src/weederTFBS.c:1050:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 3 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:1050:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 4 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:1052:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 3 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:1052:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 4 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:1054:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 3 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:1054:5: 警告: 書式 ‘%d’ は引数の型が ‘int’ であると予期されますが、第 4 引数の型は ‘long int’ です [-Wformat]
src/weederTFBS.c:458:10: 警告: warn_unused_result 属性付きで宣言されている ‘fscanf’ の戻り値を無視しています [-Wunused-result]
gcc src/weederlauncher.c -o weederlauncher.out -lm
gcc src/adviser.c -o adviser.out -O2 -lm
src/adviser.c: 関数 ‘main’ 内:
src/adviser.c:289:10: 警告: warn_unused_result 属性付きで宣言されている ‘fscanf’ の戻り値を無視しています [-Wunused-result]
gcc src/locator.c -o locator.out -O2 -lm
src/locator.c: 関数 ‘main’ 内:
src/locator.c:329:10: 警告: warn_unused_result 属性付きで宣言されている ‘fscanf’ の戻り値を無視しています [-Wunused-result]
⇒よくわからない警告文が出ていましたが、インストールはうまくいったようです。
〈Frequency fileの作成〉
Weederの場合、バックグラウンドとなるデータをあらかじめ作っておく必要があります。（転写因子の結合サイトを探す場合は、すでに用意されているファイルを使えばよいみたいです。）
（１）下記のサイトへもう一度行き、左側の「Additional tools」のFrequency maker (for adding new species to Weeder)をクリック。
http://159.149.160.51/modtools/
（２）「How to build a new frequency file for Weeder」というページに移動するので、上から２つ目の「Download this program (works under UNIX/LINUX/MacOSX/Cygwin)」を右クリック、リンク先のファイルを保存してください。（名前は「wfrequency_maker.cpp」のままでダウンロードする。）
（３）ファイルをダウンロードしたディレクトリに移動し、コンパイルを実行。
$ g++ wfrequency_maker.cpp -o wfrequency_maker
（４）バックグラウンドとして使用する配列のFASTAファイルのパス、およびバックグラウンドデータを呼び出す際に用いる任意の２文字を入力して、Frequency fileを作成する。
$ ./wfrequency_maker ./任意のパス/INPUT_file.fasta XX（任意の２文字）
Reading sequences file (avoiding duplicates)... done - Read sequence(s): 9360 Good Sequences: 9325 Duplicated sequence(s): 34
Building 6mers list... done
Finding exact occurrences... done
Counting 6mers occurrences...4096/4096 done
Writing freq file: UP.6.freq done
Building 8mers list... done
Finding exact occurrences... done
Counting 8mers occurrences...65527/65527 done
Writing freq file: UP.8.freq done
XX.6.freq
XX.8.freq
⇒上記の２種類のファイルができるので、これらのファイルをWeeder1.4.2ディレクトリ中の「FreqFiles」へ移動させる。
〈モチーフ検索〉
$ ./weederlauncher.out ./任意のパス/INPUT_file.fasta XX（バックグラウンドデータ） analysis_mode
⇒上記はあくまでコマンドの入力例になります。詳しくは、「weedermanual.pdf」を呼んでみてください。
・パラメータ
analysis_modeはいくつか用意されています。
small: 6mer、8merのモチーフを検索する
median: 6mer、8mer、10merのモチーフを検索する
large: 6mer、8mer、10mer、12merのモチーフを検索する
などなど。
・オプションパラメータ
M: モチーフ配列が同一の配列に対して「２回以上」出現するモチーフ配列を検索する。
T[1-n]: モチーフ検索の回数を指定する。（1-n回: default 10）
などなど。
〈参考文献〉
Weeder Web: discovery of transcription factor binding sites in a set of sequences from co-regulated genes
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC441603/
Motif Discovery （BioProspectorを使ってみた）
ChIP-seqやCLIP-seqを行うことにより、転写因子が結合する領域や、RNA結合タンパク質が結合する領域を絞り込むことができます。
さらに、それらの結合領域の配列データから、転写因子やRNA結合タンパク質が結合するであろうコンセンサス配列（モチーフ配列）を調べるというのが一般的な解析の流れだと思います。
ChIP-seqやCLIP-seq用に開発された解析ツールも存在しますが、今回は古典的なモチーフ検索を行うツールによるコンセンサス配列（モチーフ配列）の予測を行なってみたいと思います。
これらのモチーフ検索ソフトを利用することによって、ChIP-seqやCLIP-seqの解析のみならず、3'UTRやintron中などにひそむコンセンサス配列（モチーフ配列）を予測することができます。（予測精度は定かではありませんが…。）
〈インストール方法〉
（１）BioProspector 2004 releaseをダウンロードする。
（２）ダウンロードしてきたファイルを解凍する。
$ unzip Bioprospector.zip
（３）Bioprospectorのディレクトリに移動する。
$ cd Bioprospector
（４）LinuxをOSとして使用している場合は、「BioProspector.linux」のファイル名を「BioProspector」に変更する。
$ mv BioProspector.linux BioProspector
さらに、それらの結合領域の配列データから、転写因子やRNA結合タンパク質が結合するであろうコンセンサス配列（モチーフ配列）を調べるというのが一般的な解析の流れだと思います。
ChIP-seqやCLIP-seq用に開発された解析ツールも存在しますが、今回は古典的なモチーフ検索を行うツールによるコンセンサス配列（モチーフ配列）の予測を行なってみたいと思います。
これらのモチーフ検索ソフトを利用することによって、ChIP-seqやCLIP-seqの解析のみならず、3'UTRやintron中などにひそむコンセンサス配列（モチーフ配列）を予測することができます。（予測精度は定かではありませんが…。）
〈インストール方法〉
（１）BioProspector 2004 releaseをダウンロードする。
（２）ダウンロードしてきたファイルを解凍する。
$ unzip Bioprospector.zip
（３）Bioprospectorのディレクトリに移動する。
$ cd Bioprospector
（４）LinuxをOSとして使用している場合は、「BioProspector.linux」のファイル名を「BioProspector」に変更する。
$ mv BioProspector.linux BioProspector
（５）パーミッションの変更。
$ chmod 555 BioProspector
〈モチーフ検索〉
$ ./BioProspector -i INPUT_file.fasta -b Background_file.fasta -W motif_width -n number_of_times_trying_to_find_motif -d 1 -o Result_file.txt
⇒上記はあくまでコマンドの入力例になります。詳しくは、「BioProspector.README」を呼んでみてください。
・オプションパラメータ
-W [4-50]: モチーフ配列の長さを決める。（4-50mer: Default 10mer）
-o Result_file.txt: 結果の出力先を指定。
-b Background_file.fasta: バックグラウンドの配列データ（ゲノム配列、解析対象となる全遺伝子の配列データ（Promoter領域、Intron、3'UTR）など。）
-n [1-200]: モチーフ検索を行う回数。（1-200回: Default 40）
-r [1-n]: スコアの高いモチーフ配列を出力する数（１以上: Default 5）
-d 1: 投げた配列のAntisense鎖は解析に使用しない。（Default 0: 使用する）
などなど。
あとは試行錯誤で、自分の解析目的に応じてパラメータをいじくる。
〈参考文献〉・オプションパラメータ
-W [4-50]: モチーフ配列の長さを決める。（4-50mer: Default 10mer）
-o Result_file.txt: 結果の出力先を指定。
-b Background_file.fasta: バックグラウンドの配列データ（ゲノム配列、解析対象となる全遺伝子の配列データ（Promoter領域、Intron、3'UTR）など。）
-n [1-200]: モチーフ検索を行う回数。（1-200回: Default 40）
-r [1-n]: スコアの高いモチーフ配列を出力する数（１以上: Default 5）
-d 1: 投げた配列のAntisense鎖は解析に使用しない。（Default 0: 使用する）
などなど。
あとは試行錯誤で、自分の解析目的に応じてパラメータをいじくる。
BioProspector: discovering conserved DNA motifs in upstream regulatory regions of co-expressed genes
2013年3月28日木曜日
ディスク容量の確認（Linux）
Linux OSを使用していると、「ディスクの容量の空きってどのくらいなんだろう？」って気になるときがあります。
そんな時、以下のコマンドを打つだけでディスク容量を把握することができます。
$df -h
Filesystem Size Used Avail Use% Mounted on
/dev/sda1 19G 5.5G 13G 31% /
udev 494M 4.0K 494M 1% /dev
tmpfs 201M 1.1M 200M 1% /run
none 5.0M 0 5.0M 0% /run/lock
none 501M 124K 501M 1% /run/shm
こんな感じ。
そんな時、以下のコマンドを打つだけでディスク容量を把握することができます。
$df -h
Filesystem Size Used Avail Use% Mounted on
/dev/sda1 19G 5.5G 13G 31% /
udev 494M 4.0K 494M 1% /dev
tmpfs 201M 1.1M 200M 1% /run
none 5.0M 0 5.0M 0% /run/lock
none 501M 124K 501M 1% /run/shm
こんな感じ。
CPANを利用してPerlモジュールをインストールする（Linux, Ubuntu）
Perlモジュールをインストールしたいとき、CPANのサイトからわざわざモジュールを一個一個ダウンロードしてインストールするのは億劫です。
また、特定のPerlモジュールをインストールする際、別のPerlモジュールも必要となる場合もあるため、必要なモジュールを一括してインストールできれば便利です。
そこで、コマンドを使った方法で、Perlモジュールをインストールする方法を紹介したいと思います。この方法を用いると、足りないPerlモジュールを追加で自動的にインストールしてくれます。
〈CPAN経由でのインストール〉
（１）rootユーザーになり、管理者権限を獲得する。
$ su
（２）CPANを呼び出す。
$ perl -MCPAN -e shell
⇒初めて場合、いろいろな質問が出てくるが、「Enterキー」を押すか「y」を入力して先に進む。
（３）Perlモジュールをインストールする。
cpan[1]> install モジュール名
（４）CPANを終了させる。
cpan[2]> q
（５）rootユーザーから抜ける。
$ exit
〈指定したPerlモジュールがすでにインストールされている場合〉
（１）インストールを実行する。
cpan[3]> install HTTP::Request
HTTP::Request is up to date (6.00).
cpan[4]> install HTTP::Request::Common
HTTP::Request::Common is up to date (6.04).
上記のように表示されます。
〈指定したPerlモジュールがうまくインストールできない場合〉
（１）インストールを実行する。
cpan[2]> install XML::Parser::Expat
Running install for module 'XML::Parser::Expat'
Running make for T/TO/TODDR/XML-Parser-2.41.tar.gz
Has already been unwrapped into directory /root/.cpan/build/XML-Parser-2.41-fxgeS3
Could not make: Unknown error
Running make test
Can't test without successful make
Running make install
Make had returned bad status, install seems impossible
（２）Perlモジュールが保存されているディレクトリに移動して原因を探る。
cpan[3]> q
$ cd /root/.cpan/build/XML-Parser-2.41-fxgeS3
$ perl Makefile.PL
$ make
…
cc -c -D_REENTRANT -D_GNU_SOURCE -DDEBIAN -fno-strict-aliasing -pipe -fstack-protector -I/usr/local/include -D_LARGEFILE_SOURCE -D_FILE_OFFSET_BITS=64 -O2 -g -DVERSION=\"2.41\" -DXS_VERSION=\"2.41\" -fPIC "-I/usr/lib/perl/5.14/CORE" Expat.c
Expat.xs:12:19: 致命的エラー: expat.h: そのようなファイルやディレクトリはありません
コンパイルを停止しました。
make[1]: *** [Expat.o] エラー 1
make[1]: ディレクトリ `/root/.cpan/build/XML-Parser-2.41-fxgeS3/Expat' から出ます
make: *** [subdirs] エラー 2
⇒どうやら、expat.hがLinuxに入っていないために、エラーが起こっているようです。
（３）expat.hをapt-getコマンドでインストール
$ apt-file search /usr/include/expat.h
libexpat1-dev: /usr/include/expat.h
$ sudo apt-get install libexpat1-dev
[sudo] password for imamachi:
（４）再度、CPANでインストールを試みる。
$ perl -MCPAN -e shell
cpan[1]> install モジュール名
また、特定のPerlモジュールをインストールする際、別のPerlモジュールも必要となる場合もあるため、必要なモジュールを一括してインストールできれば便利です。
そこで、コマンドを使った方法で、Perlモジュールをインストールする方法を紹介したいと思います。この方法を用いると、足りないPerlモジュールを追加で自動的にインストールしてくれます。
〈CPAN経由でのインストール〉
（１）rootユーザーになり、管理者権限を獲得する。
$ su
（２）CPANを呼び出す。
$ perl -MCPAN -e shell
⇒初めて場合、いろいろな質問が出てくるが、「Enterキー」を押すか「y」を入力して先に進む。
（３）Perlモジュールをインストールする。
cpan[1]> install モジュール名
（４）CPANを終了させる。
cpan[2]> q
（５）rootユーザーから抜ける。
$ exit
〈指定したPerlモジュールがすでにインストールされている場合〉
（１）インストールを実行する。
cpan[3]> install HTTP::Request
HTTP::Request is up to date (6.00).
cpan[4]> install HTTP::Request::Common
HTTP::Request::Common is up to date (6.04).
上記のように表示されます。
〈指定したPerlモジュールがうまくインストールできない場合〉
（１）インストールを実行する。
cpan[2]> install XML::Parser::Expat
Running install for module 'XML::Parser::Expat'
Running make for T/TO/TODDR/XML-Parser-2.41.tar.gz
Has already been unwrapped into directory /root/.cpan/build/XML-Parser-2.41-fxgeS3
Could not make: Unknown error
Running make test
Can't test without successful make
Running make install
Make had returned bad status, install seems impossible
（２）Perlモジュールが保存されているディレクトリに移動して原因を探る。
cpan[3]> q
$ cd /root/.cpan/build/XML-Parser-2.41-fxgeS3
$ perl Makefile.PL
$ make
…
cc -c -D_REENTRANT -D_GNU_SOURCE -DDEBIAN -fno-strict-aliasing -pipe -fstack-protector -I/usr/local/include -D_LARGEFILE_SOURCE -D_FILE_OFFSET_BITS=64 -O2 -g -DVERSION=\"2.41\" -DXS_VERSION=\"2.41\" -fPIC "-I/usr/lib/perl/5.14/CORE" Expat.c
Expat.xs:12:19: 致命的エラー: expat.h: そのようなファイルやディレクトリはありません
コンパイルを停止しました。
make[1]: *** [Expat.o] エラー 1
make[1]: ディレクトリ `/root/.cpan/build/XML-Parser-2.41-fxgeS3/Expat' から出ます
make: *** [subdirs] エラー 2
⇒どうやら、expat.hがLinuxに入っていないために、エラーが起こっているようです。
（３）expat.hをapt-getコマンドでインストール
$ apt-file search /usr/include/expat.h
libexpat1-dev: /usr/include/expat.h
$ sudo apt-get install libexpat1-dev
[sudo] password for imamachi:
（４）再度、CPANでインストールを試みる。
$ perl -MCPAN -e shell
cpan[1]> install モジュール名
2013年3月23日土曜日
Sequence Read Archive (SRA)ファイルをfastqファイルに変換
論文中で使用したRNA-seqなどのNGSデータは、NCBIのホームページの「Gene Expression Omnibus (GEO)」にすべて登録されています。
しかし、このGEOに登録されているNGSデータは見慣れた「fastq」ファイルではなく、「Sequence Read Archive (SRA)」と呼ばれるファイルとして管理されています。
そのため、これらNGSデータを使用する際には、「sra」ファイルから「fastq」ファイルへ変換する必要があります。
ここでは、「sra→fastq」への変換方法の一例を紹介したいと思います。
変換ツールは、NCBIによって配布されている「SRA toolkit」を使用します。
〈SRA toolkitのダウンロード〉
（１）NCBIのサイトの上部にある検索欄のプルダウンメニューから「SRA」を選択し、検索欄には何も入れずに「Search」をクリック。
（２）Toolsの項目中にある「SRA software」をクリック。
（３）「1. NCBI SRA Toolkit latest release (February 20 2013, version 2.3.1 release) compiled binaries and md5 checksums*:」以下から、自分の使用しているOSに対応したSRA toolkitを選択し、ダウンロード。
〈SRA toolkitのインストール〉
（１）ダウンロードしたファイルを保存したディレクトリに移動。
$ cd Download/
$ ls
sratoolkit.2.3.1-ubuntu32.tar.gz
（２）ダウンロードしたファイルを解凍する。
$ tar zxvf sratoolkit.2.3.1-ubuntu32.tar.gz
（３）コンパイルされたプログラム（コマンド）に対して、パスを通す。
⇒SRA toolkitでは、すでにコンパイル済みにプログラムが「sratoolkit2.3.1-ubuntu32/bin」中に入っているので、そのディレクトリにパスを通すだけでOKです。
$ cd sratoolkit.2.3.1-ubuntu32/bin/ #パスを通したいファイルの所在を確認
$ ls
abi-dump ncbi sra-stat
abi-dump.2 nenctool sra-stat.2
abi-dump.2.3.1 nenctool.2 sra-stat.2.3.1
abi-load nenctool.2.3.1 srf-load
abi-load.2 pacbio-load srf-load.2
abi-load.2.3.1 pacbio-load.2 srf-load.2.3.1
align-info pacbio-load.2.3.1 test-sra
align-info.2 prefetch test-sra.2
align-info.2.3.1 prefetch.2 test-sra.2.3.1
bam-load prefetch.2.3.1 vdb-config
bam-load.2 rcexplain vdb-config.2
bam-load.2.3.1 rcexplain.2 vdb-config.2.3.1
cg-load rcexplain.2.3.1 vdb-copy
cg-load.2 refseq-load vdb-copy.2
cg-load.2.3.1 refseq-load.2 vdb-copy.2.3.1
configuration-assistant.perl refseq-load.2.3.1 vdb-decrypt
fastq-dump sam-dump vdb-decrypt.2
fastq-dump.2 sam-dump.2 vdb-decrypt.2.3.1
fastq-dump.2.3.1 sam-dump.2.3.1 vdb-dump
fastq-load sff-dump vdb-dump.2
fastq-load.2 sff-dump.2 vdb-dump.2.3.1
fastq-load.2.3.1 sff-dump.2.3.1 vdb-encrypt
helicos-load sff-load vdb-encrypt.2
helicos-load.2 sff-load.2 vdb-encrypt.2.3.1
helicos-load.2.3.1 sff-load.2.3.1 vdb-lock
illumina-dump sra-dbcc vdb-lock.2
illumina-dump.2 sra-dbcc.2 vdb-lock.2.3.1
illumina-dump.2.3.1 sra-dbcc.2.3.1 vdb-passwd
illumina-load sra-kar vdb-passwd.2
illumina-load.2 sra-kar.2 vdb-passwd.2.3.1
illumina-load.2.3.1 sra-kar.2.3.1 vdb-unlock
kar sra-pileup vdb-unlock.2
kar.2 sra-pileup.2 vdb-unlock.2.3.1
kar.2.3.1 sra-pileup.2.3.1 vdb-validate
kdbmeta sra-sort vdb-validate.2
kdbmeta.2 sra-sort.2 vdb-validate.2.3.1
kdbmeta.2.3.1 sra-sort.2.3.1
$ vi ~/.profile #viエディタを用いて、パスを通すディレクトリの所在を記述
赤文字で書いた部分を追加。（パスは任意にディレクトリに対して。）
export PATH="$PATH:/home/imamachi/file/samtools-0.1.19:/home/imamachi/Download/sratoolkit.2.3.1-ubuntu32/bin/"
$ source ~/.profile #上記の変更を適応
〈sra→fastqへの変換〉
（１）変換したいsraファイルのディレクトリへ移動後、下記のcommandを実行し、fastqファイルへの変換。
$ fastq-dump ./SRR317197.sra -O ./sra_result
⇒コマンドライン引数: 変換したいsraファイルの名前/ディレクトリの指定
⇒-O: 変換したfastqファイルの保存先のディレクトリの指定
以上。
しかし、このGEOに登録されているNGSデータは見慣れた「fastq」ファイルではなく、「Sequence Read Archive (SRA)」と呼ばれるファイルとして管理されています。
そのため、これらNGSデータを使用する際には、「sra」ファイルから「fastq」ファイルへ変換する必要があります。
ここでは、「sra→fastq」への変換方法の一例を紹介したいと思います。
変換ツールは、NCBIによって配布されている「SRA toolkit」を使用します。
〈SRA toolkitのダウンロード〉
（１）NCBIのサイトの上部にある検索欄のプルダウンメニューから「SRA」を選択し、検索欄には何も入れずに「Search」をクリック。
（２）Toolsの項目中にある「SRA software」をクリック。
（３）「1. NCBI SRA Toolkit latest release (February 20 2013, version 2.3.1 release) compiled binaries and md5 checksums*:」以下から、自分の使用しているOSに対応したSRA toolkitを選択し、ダウンロード。
〈SRA toolkitのインストール〉
（１）ダウンロードしたファイルを保存したディレクトリに移動。
$ cd Download/
$ ls
sratoolkit.2.3.1-ubuntu32.tar.gz
（２）ダウンロードしたファイルを解凍する。
$ tar zxvf sratoolkit.2.3.1-ubuntu32.tar.gz
（３）コンパイルされたプログラム（コマンド）に対して、パスを通す。
⇒SRA toolkitでは、すでにコンパイル済みにプログラムが「sratoolkit2.3.1-ubuntu32/bin」中に入っているので、そのディレクトリにパスを通すだけでOKです。
$ cd sratoolkit.2.3.1-ubuntu32/bin/ #パスを通したいファイルの所在を確認
$ ls
abi-dump ncbi sra-stat
abi-dump.2 nenctool sra-stat.2
abi-dump.2.3.1 nenctool.2 sra-stat.2.3.1
abi-load nenctool.2.3.1 srf-load
abi-load.2 pacbio-load srf-load.2
abi-load.2.3.1 pacbio-load.2 srf-load.2.3.1
align-info pacbio-load.2.3.1 test-sra
align-info.2 prefetch test-sra.2
align-info.2.3.1 prefetch.2 test-sra.2.3.1
bam-load prefetch.2.3.1 vdb-config
bam-load.2 rcexplain vdb-config.2
bam-load.2.3.1 rcexplain.2 vdb-config.2.3.1
cg-load rcexplain.2.3.1 vdb-copy
cg-load.2 refseq-load vdb-copy.2
cg-load.2.3.1 refseq-load.2 vdb-copy.2.3.1
configuration-assistant.perl refseq-load.2.3.1 vdb-decrypt
fastq-dump sam-dump vdb-decrypt.2
fastq-dump.2 sam-dump.2 vdb-decrypt.2.3.1
fastq-dump.2.3.1 sam-dump.2.3.1 vdb-dump
fastq-load sff-dump vdb-dump.2
fastq-load.2 sff-dump.2 vdb-dump.2.3.1
fastq-load.2.3.1 sff-dump.2.3.1 vdb-encrypt
helicos-load sff-load vdb-encrypt.2
helicos-load.2 sff-load.2 vdb-encrypt.2.3.1
helicos-load.2.3.1 sff-load.2.3.1 vdb-lock
illumina-dump sra-dbcc vdb-lock.2
illumina-dump.2 sra-dbcc.2 vdb-lock.2.3.1
illumina-dump.2.3.1 sra-dbcc.2.3.1 vdb-passwd
illumina-load sra-kar vdb-passwd.2
illumina-load.2 sra-kar.2 vdb-passwd.2.3.1
illumina-load.2.3.1 sra-kar.2.3.1 vdb-unlock
kar sra-pileup vdb-unlock.2
kar.2 sra-pileup.2 vdb-unlock.2.3.1
kar.2.3.1 sra-pileup.2.3.1 vdb-validate
kdbmeta sra-sort vdb-validate.2
kdbmeta.2 sra-sort.2 vdb-validate.2.3.1
kdbmeta.2.3.1 sra-sort.2.3.1
$ vi ~/.profile #viエディタを用いて、パスを通すディレクトリの所在を記述
赤文字で書いた部分を追加。（パスは任意にディレクトリに対して。）
export PATH="$PATH:/home/imamachi/file/samtools-0.1.19:/home/imamachi/Download/sratoolkit.2.3.1-ubuntu32/bin/"
$ source ~/.profile #上記の変更を適応
〈sra→fastqへの変換〉
（１）変換したいsraファイルのディレクトリへ移動後、下記のcommandを実行し、fastqファイルへの変換。
$ fastq-dump ./SRR317197.sra -O ./sra_result
⇒コマンドライン引数: 変換したいsraファイルの名前/ディレクトリの指定
⇒-O: 変換したfastqファイルの保存先のディレクトリの指定
以上。
Synapticパッケージマネジャー（Linux OS, Ubuntu）
Ubuntu（Linux OSの１つ）でのみ利用できる便利ツールの紹介です。
Ubuntuには「Synapticパッケージマネジャー」と呼ばれるグラフィカルなパッケージ管理ツールをインストールすることができます。
このツールを使うことで、ソフトウェアのダウンロード・インストールから、ソフトウェアのアップデート、削除まで行うことができます。（もちろん、パスも自動的に通してくれます。）
バイオ系のソフトウェアもこのSynapticパッケージマネジャーでインストール可能なので、かなり便利です。
ただ、このツール上でインストールするソフトが最新のものでない可能性があるので、注意が必要です。
〈バイオ系のソフトウェア一覧〉
・bedtools
・fastx-toolkit
・Picard
・bowtie
・cufflinks
・ncbi-blast+
などなど
※いずれもSynapticパッケージマネジャーでインストール可能なもののみ。
〈インストール方法〉
（１）Ubuntuの「Ubuntuソフトウェアセンター」で「Synapticパッケージマネジャー」で検索。
（２）候補が出てくるので、「Synapticパッケージマネジャー」を選択し、「インストール」ボタンをクリックする。
〈ソフトウェアのインストール方法〉
（１）Synapticパッケージマネジャーを起動すると、「Authentication is required to run the Synaptic Package Manager」というウインドウが最初に開くので、認証パスワードを入力。
（２）「クイック検索」欄でインストールしたいソフトの名前を入力し検索。
（３）インストールしたいにチェックをつける。
（４）上部のリボンにある「適用」ボタンをクリックするだけでインストールが自動で行われる。
Ubuntuには「Synapticパッケージマネジャー」と呼ばれるグラフィカルなパッケージ管理ツールをインストールすることができます。
このツールを使うことで、ソフトウェアのダウンロード・インストールから、ソフトウェアのアップデート、削除まで行うことができます。（もちろん、パスも自動的に通してくれます。）
バイオ系のソフトウェアもこのSynapticパッケージマネジャーでインストール可能なので、かなり便利です。
ただ、このツール上でインストールするソフトが最新のものでない可能性があるので、注意が必要です。
〈バイオ系のソフトウェア一覧〉
・bedtools
・fastx-toolkit
・Picard
・bowtie
・cufflinks
・ncbi-blast+
などなど
※いずれもSynapticパッケージマネジャーでインストール可能なもののみ。
〈インストール方法〉
（１）Ubuntuの「Ubuntuソフトウェアセンター」で「Synapticパッケージマネジャー」で検索。
（２）候補が出てくるので、「Synapticパッケージマネジャー」を選択し、「インストール」ボタンをクリックする。
〈ソフトウェアのインストール方法〉
（１）Synapticパッケージマネジャーを起動すると、「Authentication is required to run the Synaptic Package Manager」というウインドウが最初に開くので、認証パスワードを入力。
（２）「クイック検索」欄でインストールしたいソフトの名前を入力し検索。
（３）インストールしたいにチェックをつける。
（４）上部のリボンにある「適用」ボタンをクリックするだけでインストールが自動で行われる。
ファイルの解凍方法（Linux OS）
各種圧縮ファイルの解凍の仕方についてのメモ。
（１）.tar.gz (.tgz)
$ tar zxvf INPUT_file.tar.gz
（２）.tar.bz2 (.tbz)
$ tar jxvf INPUT_file.tar.bz2
（３）.gz
$ gunzip -d INPUT_file.gz
（４）.bz2
$ bzip2 -d INPUT_file.bz2
（５）.zip
$ unzip INPUT_file.zip
以上。
（１）.tar.gz (.tgz)
$ tar zxvf INPUT_file.tar.gz
（２）.tar.bz2 (.tbz)
$ tar jxvf INPUT_file.tar.bz2
（３）.gz
$ gunzip -d INPUT_file.gz
（４）.bz2
$ bzip2 -d INPUT_file.bz2
（５）.zip
$ unzip INPUT_file.zip
以上。
Pythonのインストール方法（Linux OS）
Pythonのインストール方法のメモ。
（１）Pythonのホームページに移動。
（２）左にある「DOWNROAD」をクリック。
（３）Python 2.7.3 bzipped source tarball (for Linux, Unix or Mac OS X, more compressed)を選択し、ファイルをダウンロードする。
（４）ダウンロードしたファイルの解凍
$ tar xjvf Python-2.7.3.tar.bz2
（５）解凍したディレクトリの場所に移動。
$ cd Python-2.7.3/
（６）以下のコマンドに順番に実行。
$ ./configure
$ make
$ sudo make install
（７）Pythonがインストールされているか確認。
$ python -h
usage: python [option] ... [-c cmd | -m mod | file | -] [arg] ...
Options and arguments (and corresponding environment variables):
-B : don't write .py[co] files on import; also PYTHONDONTWRITEBYTECODE=x
-c cmd : program passed in as string (terminates option list)
-d : debug output from parser; also PYTHONDEBUG=x
-E : ignore PYTHON* environment variables (such as PYTHONPATH)
-h : print this help message and exit (also --help)
-i : inspect interactively after running script; forces a prompt even
if stdin does not appear to be a terminal; also PYTHONINSPECT=x
-m mod : run library module as a script (terminates option list)
-O : optimize generated bytecode slightly; also PYTHONOPTIMIZE=x
-OO : remove doc-strings in addition to the -O optimizations
-R : use a pseudo-random salt to make hash() values of various types be
unpredictable between separate invocations of the interpreter, as
a defense against denial-of-service attacks
-Q arg : division options: -Qold (default), -Qwarn, -Qwarnall, -Qnew
-s : don't add user site directory to sys.path; also PYTHONNOUSERSITE
-S : don't imply 'import site' on initialization
-t : issue warnings about inconsistent tab usage (-tt: issue errors)
-u : unbuffered binary stdout and stderr; also PYTHONUNBUFFERED=x
see man page for details on internal buffering relating to '-u'
-v : verbose (trace import statements); also PYTHONVERBOSE=x
can be supplied multiple times to increase verbosity
-V : print the Python version number and exit (also --version)
-W arg : warning control; arg is action:message:category:module:lineno
also PYTHONWARNINGS=arg
-x : skip first line of source, allowing use of non-Unix forms of #!cmd
-3 : warn about Python 3.x incompatibilities that 2to3 cannot trivially fix
file : program read from script file
- : program read from stdin (default; interactive mode if a tty)
arg ...: arguments passed to program in sys.argv[1:]
Other environment variables:
PYTHONSTARTUP: file executed on interactive startup (no default)
PYTHONPATH : ':'-separated list of directories prefixed to the
default module search path. The result is sys.path.
PYTHONHOME : alternate <prefix> directory (or <prefix>:<exec_prefix>).
The default module search path uses <prefix>/pythonX.X.
PYTHONCASEOK : ignore case in 'import' statements (Windows).
PYTHONIOENCODING: Encoding[:errors] used for stdin/stdout/stderr.
PYTHONHASHSEED: if this variable is set to 'random', the effect is the same
as specifying the -R option: a random value is used to seed the hashes of
str, bytes and datetime objects. It can also be set to an integer
in the range [0,4294967295] to get hash values with a predictable seed.
$ python --version
Python 2.7.3
以上。
（１）Pythonのホームページに移動。
（２）左にある「DOWNROAD」をクリック。
（３）Python 2.7.3 bzipped source tarball (for Linux, Unix or Mac OS X, more compressed)を選択し、ファイルをダウンロードする。
（４）ダウンロードしたファイルの解凍
$ tar xjvf Python-2.7.3.tar.bz2
（５）解凍したディレクトリの場所に移動。
$ cd Python-2.7.3/
（６）以下のコマンドに順番に実行。
$ ./configure
$ make
$ sudo make install
（７）Pythonがインストールされているか確認。
$ python -h
usage: python [option] ... [-c cmd | -m mod | file | -] [arg] ...
Options and arguments (and corresponding environment variables):
-B : don't write .py[co] files on import; also PYTHONDONTWRITEBYTECODE=x
-c cmd : program passed in as string (terminates option list)
-d : debug output from parser; also PYTHONDEBUG=x
-E : ignore PYTHON* environment variables (such as PYTHONPATH)
-h : print this help message and exit (also --help)
-i : inspect interactively after running script; forces a prompt even
if stdin does not appear to be a terminal; also PYTHONINSPECT=x
-m mod : run library module as a script (terminates option list)
-O : optimize generated bytecode slightly; also PYTHONOPTIMIZE=x
-OO : remove doc-strings in addition to the -O optimizations
-R : use a pseudo-random salt to make hash() values of various types be
unpredictable between separate invocations of the interpreter, as
a defense against denial-of-service attacks
-Q arg : division options: -Qold (default), -Qwarn, -Qwarnall, -Qnew
-s : don't add user site directory to sys.path; also PYTHONNOUSERSITE
-S : don't imply 'import site' on initialization
-t : issue warnings about inconsistent tab usage (-tt: issue errors)
-u : unbuffered binary stdout and stderr; also PYTHONUNBUFFERED=x
see man page for details on internal buffering relating to '-u'
-v : verbose (trace import statements); also PYTHONVERBOSE=x
can be supplied multiple times to increase verbosity
-V : print the Python version number and exit (also --version)
-W arg : warning control; arg is action:message:category:module:lineno
also PYTHONWARNINGS=arg
-x : skip first line of source, allowing use of non-Unix forms of #!cmd
-3 : warn about Python 3.x incompatibilities that 2to3 cannot trivially fix
file : program read from script file
- : program read from stdin (default; interactive mode if a tty)
arg ...: arguments passed to program in sys.argv[1:]
Other environment variables:
PYTHONSTARTUP: file executed on interactive startup (no default)
PYTHONPATH : ':'-separated list of directories prefixed to the
default module search path. The result is sys.path.
PYTHONHOME : alternate <prefix> directory (or <prefix>:<exec_prefix>).
The default module search path uses <prefix>/pythonX.X.
PYTHONCASEOK : ignore case in 'import' statements (Windows).
PYTHONIOENCODING: Encoding[:errors] used for stdin/stdout/stderr.
PYTHONHASHSEED: if this variable is set to 'random', the effect is the same
as specifying the -R option: a random value is used to seed the hashes of
str, bytes and datetime objects. It can also be set to an integer
in the range [0,4294967295] to get hash values with a predictable seed.
$ python --version
Python 2.7.3
以上。
2013年3月21日木曜日
パスの通し方（Linux OS, Ubuntu）
Ubuntu（Linux OS）で起動時に自動的に読み込まれるように、パスを通してみたいと思います。
（１）現在のパスを確認。
$ echo $PATH
/usr/lib/lightdm/lightdm:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games
（２）パスを通したいコマンドのバイナリファイルの場所を確認。
$ cd /home/imamachi/file/samtools-0.1.19 #任意のファイルの場所に移動
$pwd #ファイルの場所を確認
/home/imamachi/file/samtools-0.1.19
（２）cdコマンドであらかじめルートディレクトリに移動。
$ cd
（３）viエディタを使って、「.profile」を開く。
$ vi .profile
# ~/.profile: executed by the command interpreter for login shells.
# This file is not read by bash(1), if ~/.bash_profile or ~/.bash_login
# exists.
# see /usr/share/doc/bash/examples/startup-files for examples.
# the files are located in the bash-doc package.
# the default umask is set in /etc/profile; for setting the umask
# for ssh logins, install and configure the libpam-umask package.
#umask 022
# if running bash
if [ -n "$BASH_VERSION" ]; then
# include .bashrc if it exists
if [ -f "$HOME/.bashrc" ]; then
. "$HOME/.bashrc"
fi
fi
# set PATH so it includes user's private bin if it exists
if [ -d "$HOME/bin" ] ; then
PATH="$HOME/bin:$PATH"
fi
（４）最終行に「export PATH="$PATH:/home/imamachi/file/samtools-0.1.19"」を入力
〈注意点〉
・ダブルクオーテーション「"」をつけ忘れないこと。
・＄PATHのあとに「:」を忘れないこと。
〈viエディタ入力方法〉
・主なコマンド
:q! 保存せずに終了
:wq 保存して終了
:i 現在のカーソル位置から文字を挿入
:A 現在行の末尾に文字を挿入
:x 一文字削除
:ESC コマンドモードへ移行
などなど
〈参考サイト〉
・viの使い方/基本操作
http://www.envinfo.uee.kyoto-u.ac.jp/user/susaki/command/vi.html
〈例〉
export PATH="$PATH:/home/imamachi/Download/sratoolkit.2.3.1-ubuntu64/bin:/home/imamachi/Download/samtools-0.1.19"
とか。
（５）変更をすぐに反映させる。
$ source ~/.profile
（６）変更できているか確認
$ echo $PATH
/usr/lib/lightdm/lightdm:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/home/imamachi/file/samtools-0.1.19
（７）コマンドが実際に使用できるか確認
$ samtools #任意のコマンドを入力
2013年3月20日水曜日
SAMtoolsをインストールする（Linux OS, Ubuntu）
SAMtoolsのインストールに予想以上に手間取ったので、エラーログを含めてインストール方法を書いておきたいと思います…。
〈SAMtoolsのダウンロード〉
（１）SAMtoolsのサイトへ移動。
（２）SAMtoolsサイトの「Introduction」の下のほうにある、
SAMtools is hosted by SourceForge.net. The project page is here.
The source code releases are available from the download page.
「download page」をクリック。
（３）SorceForge.netのサイトへ移動し、
Looking for the latest version? Download samtools-0.1.19.tar.bz2 (514.5 kB)
にある「Download samtools-0.1.XX.tar.bz2」をクリック。
（４）５秒待つと、自動的にダウンロードするファイルの保存場所を指定するように言ってくる。
〈SAMtoolsのインストール〉
（１）SAMtoolsのファイルを解凍。
$ tar jxvf samtools-0.1.19.tar.bz2
（２）makeコマンドでコンパイルを試みるが、致命的なエラーが発生。
$ make
make[1]: ディレクトリ `/home/imamachi/Download/samtools-0.1.19' に入ります
make[2]: ディレクトリ `/home/imamachi/Download/samtools-0.1.19' に入ります
gcc -c -g -Wall -O2 -D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE -D_USE_KNETFILE -D_CURSES_LIB=1 -DBGZF_CACHE -I. bgzf.c -o bgzf.o
In file included from bgzf.c:32:0:
bgzf.h:33:18: 致命的エラー: zlib.h: そのようなファイルやディレクトリはありません
コンパイルを停止しました。
make[2]: *** [bgzf.o] エラー 1
make[2]: ディレクトリ `/home/imamachi/Download/samtools-0.1.19' から出ます
make[1]: *** [lib-recur] エラー 1
make[1]: ディレクトリ `/home/imamachi/Download/samtools-0.1.19' から出ます
make: *** [all-recur] エラー 1
⇒どうやら「zlib.h」が入っていないせいでエラーが起こっているよう。
（３）apt-fileコマンドで「zlib.h」を検索。
$ apt-file search zlib.h
プログラム 'apt-file' はまだインストールされていません。 次のように入力することでインストールできます:
sudo apt-get install apt-file
（４）「apt-file」がインストールされていないとのことだったので、言われたとおりにコマンドを打ち、「apt-file」のインストールを試みる。
$ sudo apt-get install apt-file
[sudo] password for imamachi: #パスワードを入力
パッケージリストを読み込んでいます... 完了
依存関係ツリーを作成しています
状態情報を読み取っています... 完了
以下の特別パッケージがインストールされます:
build-essential curl dpkg-dev fakeroot g++ g++-4.6 libalgorithm-diff-perl
libalgorithm-diff-xs-perl libalgorithm-merge-perl libapt-pkg-perl
libconfig-file-perl libdpkg-perl liblist-moreutils-perl
libregexp-assemble-perl libstdc++6-4.6-dev libtimedate-perl
提案パッケージ:
debian-keyring g++-multilib g++-4.6-multilib gcc-4.6-doc libstdc++6-4.6-dbg
libstdc++6-4.6-doc
以下のパッケージが新たにインストールされます:
apt-file build-essential curl dpkg-dev fakeroot g++ g++-4.6
libalgorithm-diff-perl libalgorithm-diff-xs-perl libalgorithm-merge-perl
libapt-pkg-perl libconfig-file-perl libdpkg-perl liblist-moreutils-perl
libregexp-assemble-perl libstdc++6-4.6-dev libtimedate-perl
アップグレード: 0 個、新規インストール: 17 個、削除: 0 個、保留: 14 個。
9,640 kB のアーカイブを取得する必要があります。
この操作後に追加で 28.7 MB のディスク容量が消費されます。
続行しますか [Y/n]? y #インストールしていいか聞かれるので「y」と入力。
（５）「apt-file」がインストールできたので、「apt-file」で必要な「zlib.h」を検索。
$ apt-file search /usr/include/zlib.h #単純にzlib.hと入力するとたくさん候補が出てきてしまうので、条件を絞っておく。
zlib1g-dev: /usr/include/zlib.h
（６）「zlib1g-dev」のインストールを試みる。
$ sudo apt-get install zlib1g-dev
[sudo] password for imamachi: #パスワードを入力
（７）「zlib.h」をインストールできたので、再度、makeコマンドを実行。
$ make
make[1]: ディレクトリ `/home/imamachi/Download/samtools-0.1.19' に入ります
make[2]: ディレクトリ `/home/imamachi/Download/samtools-0.1.19' に入ります
…
bam_tview_curses.c:5:20: 致命的エラー: curses.h: そのようなファイルやディレクトリはありません
コンパイルを停止しました。
make[1]: *** [bam_tview_curses.o] エラー 1
make[1]: ディレクトリ `/home/imamachi/Download/samtools-0.1.19' から出ます
make: *** [all-recur] エラー 1
⇒今度こそうまくいくと思いきや、またしても致命的なエラーが発生(´・ω・`)
（８）しかたないので、「curses.h」のインストールを試みる。
$ apt-file search /usr/include/curses.h
libncurses5-dev: /usr/include/curses.h
$ sudo apt-get install libncurses5-dev
[sudo] password for imamachi: #パスワードを入力
パッケージリストを読み込んでいます... 完了
依存関係ツリーを作成しています
状態情報を読み取っています... 完了
以下の特別パッケージがインストールされます:
libtinfo-dev
提案パッケージ:
ncurses-doc
以下のパッケージが新たにインストールされます:
libncurses5-dev libtinfo-dev
アップグレード: 0 個、新規インストール: 2 個、削除: 0 個、保留: 14 個。
312 kB のアーカイブを取得する必要があります。
この操作後に追加で 1,334 kB のディスク容量が消費されます。
続行しますか [Y/n]? y #インストールしていいか聞かれるので「y」と入力。
（９）「curses.h」をインストールできたので、満を持してmakeコマンドを実行( ･`ω･´)
$ make
⇒致命的なエラーも起こらず、コンパイルできたよう。
（１０）「samtools」コマンドが利用できるようになっているか確認。
$ ./samtools
Program: samtools (Tools for alignments in the SAM format)
Version: 0.1.18-r580
Usage: samtools <command> [options]
Command: view SAM<->BAM conversion
sort sort alignment file
mpileup multi-way pileup
depth compute the depth
faidx index/extract FASTA
tview text alignment viewer
index index alignment
idxstats BAM index stats (r595 or later)
fixmate fix mate information
flagstat simple stats
calmd recalculate MD/NM tags and '=' bases
merge merge sorted alignments
rmdup remove PCR duplicates
reheader replace BAM header
cat concatenate BAMs
bedcov read depth per BED region
targetcut cut fosmid regions (for fosmid pool only)
phase phase heterozygotes
bamshuf shuffle and group alignments by name
⇒インストール完了ヽ(*´∀｀)ﾉ
あとは「パスを通せば」完了です！
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のダウンロード〉
（１）Integrative Genomics Viewer(IGV)のホームページへ行く。
（２）左側にある「Downloads」をクリック。
（３）「Log In」画面に移るので、
・初めてIGVをダウンロードする場合は、
（ⅰ）「To use IGV, registration is required. Click here to register」で「Click here」をクリック。
（ⅱ）Registrationに移るので、「Name」「Email」「Organization」を記入後、「Agree」をクリック。
・Registrationを済ませている場合は、
（ⅰ）「email address」を記入後、「Login」ボタンをクリック。
（３）「Downloads」画面に移るので、「Binary Distribution」のIGV_2.2.11.zipをクリックし、IGV本体をダウンロードする。
（４）また「igvtools」のigvtools_2.2.2.zip includes .genome files (225 MB) をクリックし、Genome情報をあわせてダウンロードしておく。
〈IGVの起動方法〉
（１）Zipファイルを解凍した中に、「igv.bat」というファイルがあり、そのファイルをクリックするとIGVが起動します。
〈使用するゲノムを導入〉
デフォルトではヒトゲノムである「hg18」しか入っていないので、これ以外のゲノムを使いたい場合は自身でIGVの中にファイルをインポートする必要があります。
そこで、先ほどダウンロードしてきた「igvtools」中のファイルを使います。
「igvtools」にはさまざまなゲノム情報のファイルが入っており、このファイルをIGVの中にインポートすることで目的のゲノムを使用することができるようになります。
（１）「igvtools」を解凍する。
（２）「igvtools」→「genomes」内にさまざまな生物種のゲノム情報が入っています。（場所を確認）
（３）IGVを起動し、メニューバーの「Genomes」→「Load genome from File…」をクリック。
（４）使用したいゲノム情報をクリック（hg19.genome等のファイル）
（５）「Loading Genome…」というウインドウが出るのでしばらく待つ。
（６）ロードが終わった後、メニューバーの右端下のプルダウンメニューから、先ほどインポートしたゲノムが選択できるようになっている。
〈BAMファイルの可視化〉
IGVでBAMファイルを可視化するにはいくつか条件があります。
・「Chromosome」および「Start Position」でソートする必要あり。
・マッピングしたファイル（.bam）およびIndexファイル（bam.bai）が必要。
「Tophat」でマッピングした場合、BAMファイルはソートされた状態で生成されますが、「Bowtie」でマッピングした場合、生成されるSAMファイルはソートされていないので注意が必要です。
SAMtoolsを使って、sortとindexの生成を行います。
コマンドは以下のとおり。（TophatからBAMファイルを生成した場合は（３）から。）
（１）SAMファイルからBAMファイルへの変換
$samtools view -bS INPUT_file.sam > OUTPUT_file.bam
⇒標準出力により、BAMファイルは生成されます。
（２）BAMファイルのソート
$ samtools sort INPUT_file.bam OUTPUT_file
⇒アウトプットのファイルには自動的に「.bam」の拡張子がつくのでコマンド内で記述する必要はありません。
（３）BAMファイルのindex生成
$ samtools index INPUT_file.bam(sort済み)
⇒今度はOUTPUTファイルの名前を指定する必要はありません。上記の例では、「INPUT_file.bam.bai」というindexファイルが作られます。
〈IGVでBAMファイルを可視化〉
（１）IGVを起動。（「igv.bat」ファイルをクリック）
（２）メニューから「File」→「Load from File…」をクリック。
（３）見たいBAMファイル（.bam）をクリック。
⇒このとき、あわせてindexファイル（.bam.bai）を開く必要はない。BAMファイルと同じ場所にindexファイルがあれば問題なくBAMファイルを開くことができる。
IGVを利用する主な利点としては、
・Web上で動かす「UCSC genome browser」と比較して、動作が軽快。
⇒UCSC genome browserはサーバーを介してデータのやり取りを行うので、どうしてもタイムラグが生じる。イライラする。
・マッピング後のBAMファイルなど様々な形式のファイルを簡単に可視化できる。
⇒UCSC genome browserでは、BAMファイルの場合、サーバー上にファイルをストレージする必要があり、アップロードにも時間を要する。
・クローズドな環境でデータを見れるので、セキュリティ上安全。
⇒UCSC genome browserだとどうしてもデータを外部に提示する（アップロードする）必要があり、セキュリティ上安全とはいえないかも。
といったところでしょうか。
ではIGVのダウンロード・起動方法などについての説明に移りたいと思います。
〈IGVのダウンロード〉
（１）Integrative Genomics Viewer(IGV)のホームページへ行く。
（２）左側にある「Downloads」をクリック。
（３）「Log In」画面に移るので、
・初めてIGVをダウンロードする場合は、
（ⅰ）「To use IGV, registration is required. Click here to register」で「Click here」をクリック。
（ⅱ）Registrationに移るので、「Name」「Email」「Organization」を記入後、「Agree」をクリック。
・Registrationを済ませている場合は、
（ⅰ）「email address」を記入後、「Login」ボタンをクリック。
（３）「Downloads」画面に移るので、「Binary Distribution」のIGV_2.2.11.zipをクリックし、IGV本体をダウンロードする。
（４）また「igvtools」のigvtools_2.2.2.zip includes .genome files (225 MB) をクリックし、Genome情報をあわせてダウンロードしておく。
〈IGVの起動方法〉
（１）Zipファイルを解凍した中に、「igv.bat」というファイルがあり、そのファイルをクリックするとIGVが起動します。
〈使用するゲノムを導入〉
デフォルトではヒトゲノムである「hg18」しか入っていないので、これ以外のゲノムを使いたい場合は自身でIGVの中にファイルをインポートする必要があります。
そこで、先ほどダウンロードしてきた「igvtools」中のファイルを使います。
「igvtools」にはさまざまなゲノム情報のファイルが入っており、このファイルをIGVの中にインポートすることで目的のゲノムを使用することができるようになります。
（１）「igvtools」を解凍する。
（２）「igvtools」→「genomes」内にさまざまな生物種のゲノム情報が入っています。（場所を確認）
（３）IGVを起動し、メニューバーの「Genomes」→「Load genome from File…」をクリック。
（４）使用したいゲノム情報をクリック（hg19.genome等のファイル）
（５）「Loading Genome…」というウインドウが出るのでしばらく待つ。
（６）ロードが終わった後、メニューバーの右端下のプルダウンメニューから、先ほどインポートしたゲノムが選択できるようになっている。
〈BAMファイルの可視化〉
IGVでBAMファイルを可視化するにはいくつか条件があります。
・「Chromosome」および「Start Position」でソートする必要あり。
・マッピングしたファイル（.bam）およびIndexファイル（bam.bai）が必要。
「Tophat」でマッピングした場合、BAMファイルはソートされた状態で生成されますが、「Bowtie」でマッピングした場合、生成されるSAMファイルはソートされていないので注意が必要です。
SAMtoolsを使って、sortとindexの生成を行います。
コマンドは以下のとおり。（TophatからBAMファイルを生成した場合は（３）から。）
（１）SAMファイルからBAMファイルへの変換
$samtools view -bS INPUT_file.sam > OUTPUT_file.bam
⇒標準出力により、BAMファイルは生成されます。
（２）BAMファイルのソート
$ samtools sort INPUT_file.bam OUTPUT_file
⇒アウトプットのファイルには自動的に「.bam」の拡張子がつくのでコマンド内で記述する必要はありません。
（３）BAMファイルのindex生成
$ samtools index INPUT_file.bam(sort済み)
⇒今度はOUTPUTファイルの名前を指定する必要はありません。上記の例では、「INPUT_file.bam.bai」というindexファイルが作られます。
〈IGVでBAMファイルを可視化〉
（１）IGVを起動。（「igv.bat」ファイルをクリック）
（２）メニューから「File」→「Load from File…」をクリック。
（３）見たいBAMファイル（.bam）をクリック。
⇒このとき、あわせてindexファイル（.bam.bai）を開く必要はない。BAMファイルと同じ場所にindexファイルがあれば問題なくBAMファイルを開くことができる。
2013年3月17日日曜日
Perlモジュールのインストール方法（Linux OS）
まず、ダウンロードしたPerlモジュールを解凍します。
$ tar zxvf BioPerl-1.6.1.tar.gz
Perlモジュールのデータをダウンロードすると、「Makefile.PL」というPerlのスクリプトが入っています。
このスクリプトを利用して以下のコマンドを打ち込むことで、モジュールをインストールすることができます。
$ perl Makefile.PL
$ make
$ make install
$ tar zxvf BioPerl-1.6.1.tar.gz
Perlモジュールのデータをダウンロードすると、「Makefile.PL」というPerlのスクリプトが入っています。
このスクリプトを利用して以下のコマンドを打ち込むことで、モジュールをインストールすることができます。
$ perl Makefile.PL
$ make
$ make install
2013年3月10日日曜日
Linux基本コマンド「ls」「cd」「less」
産総研で開催された「次世代シークエンサー解析入門」の復習として、その内容をまとめていこうと思います。とはいえ、更新ペースは週一ぐらいで。
まずはLinux基本コマンドについて。
ずばり覚えるべきコマンドはたったの３つ。
「ls」 「cd」 「less」 の３つです。
これらコマンドの基本的な機能は以下のとおりになります。
（１）ls
カレントディレクトリ中のファイル・ディレクトリを表示させる。
〈具体例〉
imamachi@imamachi-virtual-machine:~$ ls
Desktop Download Music Share Video
Document Image Sample.desktop Template
（２）cd
ディレクトリを移動する。
〈具体例〉
imamachi@imamachi-virtual-machine:~$ ls
Desktop Download Music Share Video
Document Image Sample.desktop Template
imamachi@imamachi-virtual-machine:~$ cd Download
imamachi@imamachi-virtual-machine:~/Download$ ls
refGene.txt
imamachi@imamachi-virtual-machine:~$ cd
Desktop Download Music Share Video
Document Image Sample.desktop Template
⇒「cd」だけ打つと、ルートディレクトリへ戻れる。
imamachi@imamachi-virtual-machine:~$ cd ..
Desktop Download Music Share Video
Document Image Sample.desktop Template
⇒「cd ..」と打つと、ひとつ上の階層のディレクトリへ移動できる。
（３）less
ファイルの中身を表示する。
〈具体例〉
imamachi@imamachi-virtual-machine:~$ ls
Desktop Download Music Share Video
Document Image Sample.desktop Template
imamachi@imamachi-virtual-machine:~$ cd Download/
imamachi@imamachi-virtual-machine:~/Download$ ls
refGene.txt
imamachi@imamachi-virtual-machine:~/Download$ less refGene.txt
790 NR_026775 chr6 + 26924771
26991753 26991753 26991753 3
26924771,26936236,26991032, 26924962,26936301,26991753,
0 LINC00240 unk unk -1,-1,-1,
…
ここまでは基本的な機能について紹介してきましたが、「ls」 「cd」 「less」には他にもオプションをつけることができます。
（１）ls -l
「ls」に「-l」のオプションをつけることによって、さらにカレントディレクトリ中のファイル・ディレクトリの詳細について表示させることができます。
〈具体例〉
imamachi@imamachi-virtual-machine:~/Download$ ls -l
合計 14256
-rw-rw-r-- 1 imamachi imamachi 1052432 2月 3 2012 liftOver
-rw-rw-r-- 1 imamachi imamachi 1052432 2月 3 2012 liftOver(1)
-rw-rw-r-- 1 imamachi imamachi 12491136 3月 10 22:38 refGene.txt
imamachi@imamachi-virtual-machine:~/Download$
（３）less -S
「less」に「-S」のオプションをつけることによって、改行コードが入っていなくてもウインドウサイズに合わせて折り返されていた文字列を、改行コードが入っているときのみ折り返して表示してくれます。（ただし、このオプションをつけると、ウインドウに入りきらなかった部分に関しては見せてくれません。）
〈具体例〉
〈具体例〉
imamachi@imamachi-virtual-machine:~/Download$ less -S refGene.txt
790 NR_026775 chr6 + 26924771 26991753
216 NM_007188 chr7 + 150725509 150744869
1152 NM_001242928 chr14 + 74353317 74398991
1645 NM_197964 chr7 + 139025877 139031065
257 NM_197962 chr1 - 193065594 193074608
135 NM_197961 chr15 - 65737997 65810035
⇒恥ずかしながら、講習会で初めて知りました。へたに文字列を折り返されると読みにくいファイルに関しては、このオプションを利用すると便利ですね。特に、よく使うBED, GTF, SAM, FASTQファイルの中身を見たいときには有効かな。
今日はここまで。
ちなみに
・VMware(R) Player 5.0.2 build-1031769
・Ubuntu 12.04 LTS
を利用して上記のコマンドを実行してます。
2013年3月2日土曜日
Pythonモジュールのインストール方法（Linux OS）
Pythonモジュールのデータをダウンロードすると、「seyup.py」というPythonのスクリプトが入っています。
そこで、このスクリプトを利用して以下のコマンドを打ち込むことで、モジュールをインストールすることができます。
$ python setup.py install
基本的にこのコマンドを打てばインストールできるはずです。
そこで、このスクリプトを利用して以下のコマンドを打ち込むことで、モジュールをインストールすることができます。
$ python setup.py install
基本的にこのコマンドを打てばインストールできるはずです。
UCSC genome browserの提供するbinary utilitiesを利用してみる
UCSC genome browserでは、いくつか便利なソフトウェアを提供しています。
これらのソフトウェアはLinux OSなどで動作するソフトで、パスの通してあるディレクトリにぶち込むだけで動作します。（例えば、"/bin"など。）
■UCSC genome browser -the directory of binary utilities-
http://hgdownload.cse.ucsc.edu/admin/exe/
これらのソフトウェアはLinux OSなどで動作するソフトで、パスの通してあるディレクトリにぶち込むだけで動作します。（例えば、"/bin"など。）
■UCSC genome browser -the directory of binary utilities-
http://hgdownload.cse.ucsc.edu/admin/exe/
2013年1月4日金曜日
UCSC genome browserからGTF, BED, FASTAなど様々な形式のファイルをダウンロードする（２）
前回に引き続き、UCSC genome browserの「Table browser」の活用についての説明です。
前回紹介したもの以外にも使えるデータや小技があるのでそれについて少し触れておきたいと思います。
Example1: 種間保存性の高い領域のピックアップ
Phastconsのデータをもとに、ヒトゲノム上における種間保存性の高い領域（哺乳類での比較）をBEDファイルでダウンロードしたいと思います。
（１）UCSC genome browserのトップページから、左側のメニューにある「Table browser」をクリック
（２）Table browserの各プルダウンメニューで以下を選択。
clade: Mammal
genome: Human
assembly: Feb. 2009 (GRCh37/hg19)
group: Comparative Genomics
track: Conservation
table: Mammal EI (phastConsElements46waysPlacental)
region: genome
output file: BED - browser extensible data
⇒オプションとして、Galaxyにチェックを入れるとGalaxyにデータが保存されます
file type returned: plain text / gzip compressed
⇒好みによるが、圧縮ファイルとして欲しい場合はgzipを選択
（３）「get output」をクリック
（４）Output refGene as BED
・「Incude custom track header」にチェックを入れると、custom trackとしてアップロードできるようにファイルを加工してくれる。「name」「description」「visibility」などの入力欄が用意されている。
・「Create one BED record per」で「Whole Gene」を選択
（５）「get BED」をクリック
⇒Google chromeだと、ブラウザ上でファイルが展開してしまうので、Galaxyを経由してファイルはダウンロードしてきたほうがよいかも。
（６）Output results to Galaxy as RefSeq Genesで「Send query to Galaxy」をクリック
（７）右側のヒストリーでジョブを確認。灰色は処理待ち、黄色は処理中、緑色は処理終了を表す。処理が終わったら、ジョブをクリックしてフロッピーディスクの形をした「ダウンロード」ボタンをクリックし、目的のファイルをダウンロード。
Example2: BEDファイル<=>GTFファイル相互変換
UCSC genome browserに登録されていないNONCODEやGENCODEなどのデータ（BEDやGTFファイルとして各サイトで入手可）は必ずしも統一されておらず、そのままのファイル形式ではTophatなどのマッピングソフトで使用できないことがあります。
そこで、まずCustom trackとしてアップロードしたデータをTable browserを介してデータを変換することにより目的のファイル形式に変換します。（他にも良いやり方があるかもしれません。）
（１）UCSC genome browserのトップページから、左側のメニューにある「Table browser」をクリック
（２）Table browserの各プルダウンメニューで以下を選択。
clade: Mammal
genome: Human
assembly: Feb. 2009 (GRCh37/hg19)
group: Custom track
track: 任意のデータ
table: 任意のテーブル
region: genome
output file: BED - browser extensible data
⇒オプションとして、Galaxyにチェックを入れるとGalaxyにデータが保存されます
file type returned: plain text / gzip compressed
⇒好みによるが、圧縮ファイルとして欲しい場合はgzipを選択
（３）あとは上記と同様。
前回紹介したもの以外にも使えるデータや小技があるのでそれについて少し触れておきたいと思います。
Example1: 種間保存性の高い領域のピックアップ
Phastconsのデータをもとに、ヒトゲノム上における種間保存性の高い領域（哺乳類での比較）をBEDファイルでダウンロードしたいと思います。
（１）UCSC genome browserのトップページから、左側のメニューにある「Table browser」をクリック
（２）Table browserの各プルダウンメニューで以下を選択。
clade: Mammal
genome: Human
assembly: Feb. 2009 (GRCh37/hg19)
group: Comparative Genomics
track: Conservation
table: Mammal EI (phastConsElements46waysPlacental)
region: genome
output file: BED - browser extensible data
⇒オプションとして、Galaxyにチェックを入れるとGalaxyにデータが保存されます
file type returned: plain text / gzip compressed
⇒好みによるが、圧縮ファイルとして欲しい場合はgzipを選択
（３）「get output」をクリック
（４）Output refGene as BED
・「Incude custom track header」にチェックを入れると、custom trackとしてアップロードできるようにファイルを加工してくれる。「name」「description」「visibility」などの入力欄が用意されている。
・「Create one BED record per」で「Whole Gene」を選択
（５）「get BED」をクリック
⇒Google chromeだと、ブラウザ上でファイルが展開してしまうので、Galaxyを経由してファイルはダウンロードしてきたほうがよいかも。
（６）Output results to Galaxy as RefSeq Genesで「Send query to Galaxy」をクリック
（７）右側のヒストリーでジョブを確認。灰色は処理待ち、黄色は処理中、緑色は処理終了を表す。処理が終わったら、ジョブをクリックしてフロッピーディスクの形をした「ダウンロード」ボタンをクリックし、目的のファイルをダウンロード。
Example2: BEDファイル<=>GTFファイル相互変換
UCSC genome browserに登録されていないNONCODEやGENCODEなどのデータ（BEDやGTFファイルとして各サイトで入手可）は必ずしも統一されておらず、そのままのファイル形式ではTophatなどのマッピングソフトで使用できないことがあります。
そこで、まずCustom trackとしてアップロードしたデータをTable browserを介してデータを変換することにより目的のファイル形式に変換します。（他にも良いやり方があるかもしれません。）
（１）UCSC genome browserのトップページから、左側のメニューにある「Table browser」をクリック
（２）Table browserの各プルダウンメニューで以下を選択。
clade: Mammal
genome: Human
assembly: Feb. 2009 (GRCh37/hg19)
group: Custom track
track: 任意のデータ
table: 任意のテーブル
region: genome
output file: BED - browser extensible data
⇒オプションとして、Galaxyにチェックを入れるとGalaxyにデータが保存されます
file type returned: plain text / gzip compressed
⇒好みによるが、圧縮ファイルとして欲しい場合はgzipを選択
（３）あとは上記と同様。
登録:
コメント (Atom)