そこで、Linuxの改行コードに変換するcommandを紹介します。
dos2unix ファイル名でLinuxの改行コードに変換可能です。データは上書きされます。
2013年9月14日土曜日
WindowsとLinux間での改行コードの違いによる問題
改行コードの違いから、Perlなどのスクリプトで処理した時、予期せぬところに改行が入ったりすることがあります。
そこで、Linuxの改行コードに変換するcommandを紹介します。
そこで、Linuxの改行コードに変換するcommandを紹介します。
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年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年5月26日日曜日
ディスクの容量に困ったらファイルを圧縮しておこう
大規模データを処理していると、あっという間にディスクの容量を食います。そのため、どうにかしてディスクの容量を抑える必要があります。昔解析したデータを順次消していくというのが手っ取り早い方法かも知れませんが、過去に解析したデータを残しておきたいケースもあると思います。
そこで、以下の圧縮率が高い「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月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
端末内でのコピー＆ペーストを可能にする（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月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日木曜日
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日土曜日
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」というウインドウが最初に開くので、認証パスワードを入力。
（２）「クイック検索」欄でインストールしたいソフトの名前を入力し検索。
（３）インストールしたいにチェックをつける。
（４）上部のリボンにある「適用」ボタンをクリックするだけでインストールが自動で行われる。
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
⇒インストール完了ヽ(*´∀｀)ﾉ
あとは「パスを通せば」完了です！
登録:
投稿 (Atom)