ネットブックで行うNGS解析004：ネットブックで行うRNA-seq(3)

　このコーナーではcufflinksによる発現量の推定＆遺伝子との対応付けを行う。

　念のためであるが、cufflinksのインストール＆セットアップは、まずはネットブックによるRNA-seqデーター解析(1)を参考にしてみるとよい。cufflinks, BoostC++, Eigen, samtoolsのインストールが併せて必要である (*)。

　さてネットブックで行うRNA-seqデーター解析(2)で行ったように、今シェアフォルダに以下のようなファイルが出来ているはずである。


SRR886461.sort.bam
SRR886462.sort.bam
SRR886463.sort.bam
SRR886464.sort.bam

これらソートされたBAMファイルをつかい今回は各遺伝子の発現量の推定を行う。

これはcufflinksと呼ばれるソフトウエアパッケージを用いて行われ、
このうち
　
　一つのファイルから絶対的な発現量を推定するcufflinksと、
　複数のサンプルを比較するcuffdiff

が主に使われるものである。

cufflinksについてはアラフォーからのハーバード留学研究編009:ド素人(ウエット系）のcufflinks解析(1)のところで少し述べたので、

今回はcuffdiffについて中心的にのべてみる。
大概は実験群とコントロール群の比較をすることがおおいので、こちらの方が重要かもしれない。
　
　基本的には

$cuffdiff （オプション）reference.gtf sample1.bam, sample2.bam,...., sampleN.bam　control1.bam, control2.bam...controlM.bamのように行うとよい。

オプションの部分はなくても動くが、アウトプットフォルダを決めたりやサンプルの名前を入力できたり、した方がよいので、オプションの指定方法についてもマニュアルを見ておいた方がよいかもしれない。

主なものは

-p コンピューターのコア数
-o output folder　：アウトプットを保存するディレクトリ名
-L label1, label2...　：サンプルの仮名（ラベル）設定　

くらいである。

ラベルの設定がなぜかうまく設定できないので、(アラフォーからのハーバード留学研究編010:ド素人(ウエット系）のcufflinks解析(2)参照)今回はラベル設定はしないこととした（デフォルトでq1,q2,q3,q4となる）(**)

実際のオペレーションは

$ cd /mnt/hgfs/shared
sharedフォルダにうつって

$ cuffdiff -p 2 UCSC_mm9.gtf -o LT_HSC_VS_ST_HSC_cuffdiff_result SRR886461.sort.bam, SRR886462.sort.bam SRR886463.sort.bam, SRR886464.sort.bam
とうちこめばよい。基本は上気つのBAMファイルを、染色体情報と発現量を関連付けて、LT_HSC_VS_ST_HSC_cuffdiff_resultというフォルダに出力するというものである。

するとトータル5－6時間の計算が行われ(***)

Processed 21125 loci.
Performed 86516 isoform-level transcription difference tests
Performed 79102 tss-level transcription difference tests
Performed 74672 gene-level transcription difference tests
Performed 70942 CDS-level transcription difference tests
Performed 6057 splicing tests
Performed 3600 promoter preference tests
Performing 5416 relative CDS output tests
Writing isoform-level FPKM tracking
Writing TSS group-level FPKM tracking
Writing gene-level FPKM tracking
Writing CDS-level FPKM tracking


のようなログメッセージが出れば一応成功と言ってもよいのではないだろうか？

またLT_HSC_VS_ST_HSC_cuffdiff_resultフォルダには

cds_exp.diff
cds.diff
cds.fpkm_tracking
gene_exp.diff
genes.fpkm_tracking
isoform_exp.diff
isoforms.fpkm_tracking
promoters.diff
splicing.diff
tss_group_exp.diff
tss_groups.fpkm_tracking

というファイルが出来ているはずである。

テキストファイルなので、今後はこれをエクセルで読み込んで処理するか、Ｒなどを用いて処理することが可能である。ここまでくれば大体のプロセスは終わったようなものである。(ネットブックで行うRNA-seqデーター解析(3.5)に続く）


(*)もう少しだけ詳しい説明は以下を参照
アラフォーからのハーバード留学IT編006：バカチョンBoostC++ビルト
アラフォーからのハーバード留学IT編007：サルでもできるcufflinksセットアップ
も参照のこと。

(**)そのうち解決法を見つけたい

(***)寝る前にセットして、朝結果が出ている感じなので、早くはないものの、使えなくはないスピードである。正直なところ29800円のネットブックなんかでこんなことができてしまうのはかなり衝撃である。

ネットブックで行うNGS解析003：ネットブックで行うRNA-seq(2)

　ネットブックで行うRNA-seqデーター解析(1)に引き続き、ここでは、シーケンス生データーを発現量解析であるcufflinksにかけられる形（ソートされたBMAファイルに）に変換することを試みる。
　生データーに染色体情報を関連付け、さらに染色体の順番に並べる作業である。

具体的には

1)Bowtieによるマッピング（fastq->SAM変換）
2)SAM→BAM変換
3)BAMソート

の３つのステップからなる。

1)Bowtieによるマッピング（fastq->SAM変換）
このステップはTophatというソフトで行う方がより正確なようであるが、Tophatを試みたところ今のシステムではうまくいかなかった(アラフォーからのハーバード留学研究編013：何となくTophatインストール)。メモリ不足が原因というbiopapyrusさんのブログの指摘もあり、BowtieもしくはBowtie2だというまくいくのでBowtieをこころみることととする。

　Bowtieのインストールはネットブックで行うRNA-seqデーター解析(1)で述べたように、バカチョンで
$ sudo aptitude install bowtie (*)
でOKなはずである。

次に実際にbowtieを行う上で、マッピングに必要なファイルをsharedフォルダ(**)に用意することが必要である


(１）インデックス化されたゲノム情報（ebwtファイル）

　主な生物種であれば、自分で作る必要性がない。ここからダウンロードする）
もしなければ、ゲノム情報（ここからダウンロード）を使ってbowtie-buildを行う。

(２）fastaフォーマットのゲノム情報（リファレンスゲノム。ここからダウンロード　UCSC_㎜9.faなど)。

今回はUCSCの㎜9を使用することとした（Ensemble由来だと相性が悪いこともあるのでこのシリーズ（㎜9や㎜10）が一番無難。

indexファイルはこの場合㎜9.1.ebwtのような名前（インデックス名）がついているので、faファイルもそれに準じて名前を変えておいた方がよい（mm9.faなど）

(３）RNA-seqの生データーであるfastqファイル
ネットブックで行うRNA-seqデーター解析(1)で用意した

SRR886461.fastq
SRR886462.fastq
SRR886463.fastq
SRR886464.fastq

の4つのファイルである。
である。

実際のBowtieのオペレーションはbiopapyrusさんのサイトを参照するとよい。

一番簡単な形で行うと

$ bowtie オプション インデックス　FASTQファイル 出力SAMファイル
となる。

実際にはマウスの実験データーであるSRR88461.sraファイル(シングルリードである)をベースにつくったfastqファイルをマッピングしてみると

$ bowtie -S -p 2 ㎜9 SRR886461.fastq SRR886461.sam

と打ち込むと計算が始まり、30分くらいで結果がである。一応OKそうな結果である。なお-SはSAMファイルを出すためのオプション、-pはコア数である。サンプルがペアエンドの場合やSOLiDのcsfastaの場合は、biopapyrusさんのサイトを参照のこと。

またbowtieのピットフォールについて、後で気が付いたことを述べておいた。もしbowtieがスムーズにいかないときは、もしかするとこういうことかもしれない。アラフォーからのハーバード留学研究編025：Bowtieのピットフォール参考のこと。

２）SAMをBAMに変換する
本来は１）のステップで出てきたSRR886461.samをソートしてそのままcufflinksに流せるはずなのであるが、なぜかうまくいかないので(*)、ここでsamtoolsを使いSAM->BAM変換を行う。この方が情報量が少なくなるので後々便利という話もある(**）

実際には

$ samtools view -bS XXX.sam > XXX.bam
と打ち込めばよいから

$ samtools view -bS SRR884661.sam > SRR886461.bam
のように行う。

詳しくはNGS surfur's wiki参照。

これで変換が終わったのでこのBAMファイルを染色体番号順に並び替えるBAMソートを行う。

3)BAMソートも同じくsamtoolsを用いて行う。

これもNGS surfur's wikiに書かれている通り

$ samtools sort XXX.bam XXX.sort
とやればよいので

$ samtools sort SRR886461.bam SRR886461.sort
のようにとすればよい。(***）

最終的に
SRR886461.sort.bam
SRR886462.sort.bam
SRR886463.sort.bam
SRR886464.sort.bam
というファイルが４つできているはずである。

これでこのステップは完了である。
（ネットブックで行うRNA-seqデーター解析(3)に続く）

(*)本来は
sort -k 3,3 -k 4,4n SRR886461.sam > SRR886461.sorted.sam
(http://cufflinks.cbcb.umd.edu/manual.html#input)でソートできるはずでこれをもとに、cufflinksがかけられるはずであるが、biostarの記事のようになぜか

Error: this SAM file doesn't appear to be correctly sorted! current hit is at Contig10010|nucleolin-:26, last one was at Contig1000|---NA---:198 Cufflinks requires that if your file has SQ records in the SAM header that they appear in the same order as the chromosomes names in the alignments. If there are no SQ records in the header, or if the header is missing, the alignments must be sorted lexicographically by chromsome name and by position.

といったエラーメッセージがである。

(**)漫画時々研究者さん参照。BAMの方が容量も小さいし、後の解析に色々便利とある。

(***)$ samtools sort SRR886461.bam SRR886461.sort.bamと間違って打ち込まない。SRR886461.sort.bam.bamというアウトプットファイルになる。

ネットブックで行うNGS解析002：ネットブックで行うRNA-seq(1)

　バイオインフォマティックスにはド素人の、ウエットバイオロジスト＆アナログ系のアラフォーであるが、先日来ひょんなことからバイオインフォの勉強にはまってしまうようになってしまった。

　ネットで調べながら行うと意外なことに予想以上の学習効果があり日々進歩があること、また29800円のエイサーネットブックでありながら結構いろいろなことができるためだ。
　
　実際使っている環境は

PC:Acer Aspire One(AO756-H　14C/K)

CPU:Intel, Celeron CPU847@1.47GHz　(コア数2/スレッド数2)

RAM:4.00GB

HDD:300GB

OS：Windows　8(64bit)

である。

　そこでネットブックで行うRNA-seq解析について、自分の勉強方々、サルでもできるようにマニュアル化してみることにした。(専門的な知識はNGS sufer's wikiなど参照してください)

　実際に行う過程は

０）準備段階

１）NCBIサイトのSRAサイトからRNA-seqデーター（sraファイル）をダウンロード

２）sraファイルを、汎用性の高いfastqファイルに変換：sratoolkit

３）生データーの解析(1):シーケンスデーターを染色体上にマッピング：Bowtieなど

４）生データーの解析(2)：発現量の推定＆遺伝子との対応付け：samtools, cufflinksなど

５）発現解析：Gene Ontology(GO)解析、GSEA解析など

の5つからなる。

このブログ記事シリーズ「ネットブックで行うRNA-seqデーター解析」では
主に1-4の部分を、SRAサイトにあるデーターを使って行う。

０）準備段階であるが、
（１）Window+Linux環境を整える
（２）必要なソフトウエアをインストールすること
が必要になってくる。

（１）Window+Linux環境を整える
　bowtie, samtools, cufflinksなどRNA-seqデーターを扱うソフトウエアは、MacかLinuxをベースのものが多い。このためMac（コマンドラインが動かせないといけない）もしくはLinuxマシンでないといけない。

　筆者はwindowsユーザーであるため、Linuxを動かす方法をネットで調べたところ仮想マシンでLinuxをwindow上で動かすのが一番簡単なようである。このシステムではWindowシステムと仮想Linuxシステム（VMplayer+Ubuntu）のshare folderが設定できるので、両方のシステムでファイルがいじれて非常に使い勝手がよい。初心者的はこれにいろいろRNA-seqデータを出し入れして解析するのが便利である。

ちなみにはshare folderはデフォルト（というか以下の埼玉大学の後藤 祐一先生のブログの通りにやれば）は

Windows側(shareフォルダの置きかたにもよるが)からは、
C:\User\VertualMachines\shared

Linux側からは/mnt/hgfs/shared　
（これはルートディレクトリにあって、なぜかuserディレクトリにない。userディレクトリは/computer name(home)/VertualMachine name(user）になっている気がする）

に存在する。


★仮想化Linux(Windows+Linux)については以下を参照して取り組んでほしい。

　特に埼玉大学の後藤 祐一先生のブログを参照するとすぐに使えるようになる。ここではshare folderおよびコマンドaptitudeを把握するのがポイントである。筆者のネットブックではコア数2、使用メモリ２ＧＢ、使用HDD 100GBで設定し安定して動いている。

アラフォーからのハーバード留学IT編005:アナログ人間のWindows＋Linux

アラフォーからのハーバード留学IT編008:超簡単windows8+Linux

なおLinuxコマンドについては必要なものをその都度覚えていくのがよさそうである。心配な人はまずはLinux ちょー入門でも覚えておくとよい。また基本的には大文字と小文字を区別するシステムなので、それについても注意が必要である。

よくつかうものについては少しのべておく
(１）sudo
(２）aptitude
(３）cd

(１）sudoはシステムの管理者の権限で行うコマンドであり、ソフトウエアのインストールやそのセットアップに使うことがおおい。
$sydo aptitude install bowtieのように、$sudo 個々のコマンド名で使用されることが多い。またそのコマンドにはパスワードを必要とし、そのパスワードはUbuntuセットアップ時のものをつかう。

（２）aptitudeはソフトウエアのインストールをバカチョンでやってくれるコマンド（ソフトウエア名）である。非常に便利なのでUbuntuセットアップ後にインストールした方がよい。aptitudeのインストールは、
端末を起動し

$ sudo apt-get install aptitude

と打ち込めばOKである。
これがインストールされると、ほかのソフトウエアをインストールするには

$ sudo aptitude search キーワード

で探すといい。例えばBowtie入れたいときは、

$ sudo aptitude search bowtieと打ち込めばキーワードが出てくる。

またインストールも簡単で、

$ sudo aptitude install ソフトウェアパッケージ名（aptitude searchで表示された名前）

で一瞬でインストールできる。

　（３）cd
フォルダからフォルダへ（ディレクトリからディレクトリ）移動するコマンドで一番よく使うコマンドではないかと思う。
$ cd (userディレクトリ（端末で一番最初にでてくるディレクトリにかえる)
$ cd ..　（ひとつ前のディレクトリにいく）
$ cd /　（ルートディレクトリにいく）
$ cd ~/(userディレクトリ（端末で一番最初にでてくるディレクトリにかえる)
$ cd /mnt/hgfs/shared (Windowsとのシェアフォルダに行く）
を覚えておけばことは足りるであろう。

（２）必要なソフトウエアをインストールする

★Bowtieのインストール＆セットアップについては以下を参照してセットアップしておいてほしい。基本は

$ sudo aptitude install bowtie
で一発である。

アラフォーからのハーバード留学研究編014：最後の難関?bowtieをクリア


★samtools, cufflinksインストール＆セットアップについては以下のサイトを参照のこと。

BoostC++, cufflinks, samtools, eiganのインストール＆セットアップが日必要である。

基本的には以下のの4つで一発である。

$ sudo apt-get install libboost-dev
（BoostC++インストール）

$ sudo aptitude install cufflinks
（cufflinksインストール）　

$ sudo aptitude install samtools
(samtoolsインストール）　

$ sudo aptitude install libeigen3-dev　
（eiganライブラリインストール）

アラフォーからのハーバード留学IT編006：バカチョンBoostC++ビルト

アラフォーからのハーバード留学IT編007：サルでもできるcufflinksセットアップ

次に
１）NCBIサイトのSRAサイトからRNA-seqデーター（sraファイル）をダウンロードを行う。

　RNA-seqに限らず、Chip-seq、MeDIP-seqなど次世代シークエンサーのデーターは、NCBIのSequence Read Archive (SRA)サイトから入手可能である。キーワード検索して興味のあるデータを探してみるとよいであろう。

　ここでは最も未分化な血液幹細胞であるLT-HSCとやや分化した血液幹細胞であるST-HSCのRNA-seqデーターを使ってみることにする。

SRAサイトに移動。hematopoieticで検索。ヒットした

LT-HSC1:SRR886461.sra
LT-HSC2:SRR886462.sra
ST-HSC1:SRR886463.sra
ST-HSC2:SRR886464.sra

からSRR886461.sra～SRR886464.sraの４つのファイルをダウンロード(時間がかかります)。

Lhaplusか何かで解凍するとよい（windowsで行うのがはやいです）。

そして
２）sraファイルを、汎用性の高いfastqファイルに変換：sratoolkit
の過程をおこないます。

これはNCBIのサイトにおいてある生データのファイル形式がSRAファイルとなっており、解析するためにはこれを汎用性のあるfastaqファイルに変換しないといけないためである。この変換にはNCBIの提供するsra-toolkitをダウンロードする必要がある。

　Linux版でもできないことはないが、パスを通すのがちょっと面倒なので、初心者はwindows版をダウンロードして、これでsra->fastq変換を行うのがよいかもしれない。
　
　ダウンロードしたファイルをLhaplusか何かで解凍すると、sratoolkit.2.3.3-3-win64というフォルダができる。その中のbinフォルダの中に、fastq-dumpというアイコンがあるのでこれにsraファイルをドラッグ＆ドロップすると解析が始まる（これも結構時間がかかり１ファイルあたり数時間かかる）。

　終わるとSRR886461.fastqという感じのファイルができる（5-6GBくらいか？）。これはシーケンスデーターをfasta形式でいれたものである。
これをもとに染色体への関連付け（マッピング）を行うのが次のステップである。

なおwindows版はペアエンドの場合も一つのfastqファイル(本来は2つになるべき)になってしまうため、ペアエンドのデーターの場合は使えないことに注意。　

Linux版のインストールはアラフォーからのハーバード留学研究編011：とまどいのsra -> fastq変換の脚注に書いたので参照のこと。

最後に参考になるサイトを重複も含めてのべておく

1）Linux 基本コマンド
Linux ちょー入門
Ratポータル　基本コマンド集

2）Samtoolsの使い方
NGS Sufer's Wiki

3）Cufflinkの使い方
Wolf EarsさんのTophat・Cufflinksを用いた遺伝子発現解析の方法 (4)
Ken Osakiさんの発現解析パイプラインを作るぞ！　その1：　TopHat　の使い方
Ken Osakiさんの発現解析パイプラインを作るぞ！　その2：　TopHat　の使い方 2

4)Reference.gtfの入手サイト
Broad InstituteのFTPサイト
　
（ネットブックで行うRNA-seqデーター解析(2)に続く）

ネットブックで行うNGS解析001：ショボPCでどこまでできるか？

ネットブックで行うNGS解析としたが、スパコンとかではなく家庭用のPCでもある程度はNGS解析ができる気がしてきている。

その最たるものは

１）お家でできる次世代シーケンサ解析

２）STAP細胞関係のゲノムデータを解析してみた

てな感じであるが、
みんな結構エイサーネットブックなどのショボPCで結構ゴリゴリやっているみたいだ(*)。「わらわもではやってみようぞよと乗り出した」のがこの企画である。

最近いくつか解析してみているが、できるものもあるが、やはり歯が立たないものもあると気が付き始めた。
ただ試しにトライしてみることは悪くない。これまでにも書いた記事を、ここでも再録して見やすくしてみた。

(*)
次世代シーケンス解析に関する質問掲示板Seqansers.comには、

Bowtie out of memoryの対処法を質問している人が、

Whats wrong here and what can I do without changing computer?
I am aligning onto mm9, ebwt-files are downloaded from Bowtie-website.
I have an Acer Aspire, 3GB memory. Ubuntu.

と悲痛な叫びを、述べていた。駄目な時は、駄目です。ハイ。