そこで、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ファイルの名前を指定。
登録:
投稿 (Atom)