表紙(FrontPage) | 編集(管理者用) | 差分 | 新規作成 | 一覧 | RSS | 検索 | 更新履歴

RNA-seqデータの分析について勉強する -

目次

RNA-seq は当たり前に行われるようになった。様々なデータ解析用のソフトウェア、わかりやすい日本語の解説も多数公開されている。それらのおかげで私もその成果を使える時期が到来した(もっと早く取りかかってもよかった)。 そこで、公開されているデータを、自分の研究に役立てることを目標とする。 Counter

RNA-seq データを分析する方法について解説を公開しているページ、雑誌記事、本

これらの解説を見て勉強する。

そのほか多数

どのようなデータが公開されているか

http://www.ebi.ac.uk/arrayexpress/   には、マイクロアレイの結果と RNA-seq の結果が納められている。私の研究の分野でもすでに多数のデータがある。

E-GEOD-48661  The Arabidopsis Zinc Finger Protein 3 integrates ABA and light signaling in seed germination and plant development   というのがある。これをまず分析してみたい。

https://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-48661/?keywords=arabidopsis+seed&organism=&array=&exptype[]=&exptype[]=%22sequencing+assay%22&page=1&pagesize=25

このデータは、「Comparison of ZFP3 overexpressing and wild type Arabidopsis seedlings」と書いてある。rRNA を、RiboMinus Plant Kit というのもので取り除いている。SOLID 5500xl という機械で分析している。


ArrayExpress を見てみると、イルミナ社のシーケンサーを使用している例が多い。この場合、サンプル調製には「スタンダードmRNA キット」「Total RNA キット」「Small RNA キット」の3種類がある。たいていのものはスタンダードキットを用いていて、核ゲノムから転写されポリAが付いている標準的な mRNA のみが分析される。 その場合、オルガネラゲノムから転写される mRNA などは分析できない。これはとても残念なことであるが、しかたがない。オルガネラゲノムからの転写産物を分析したければ、自分で Total RNA キットを用いてデータを取るしかないかもしれない。

しかしオルガネラゲノムから転写される mRNA の分析を実際にしている人もいるらしい。   http://bfg.oxfordjournals.org/content/12/5/454   RNA-Seq data: a goldmine for organelle research    David Roy Smith   Briefings in Functional Genomics Volume 12, Issue 5Pp. 454-456   オルガネラゲノムからの転写産物は AT-rich なことが多いので、ポリA 選択で混入することがあるらしい。 またオリゴ dT をプライマーにする場合、ポリ A がない RNA でも A リッチな部分にプライマーがついて DNA 合成が起きる可能性がある。 そういう現象を利用するので、オルガネラゲノムからの転写産物にたいしては定量的なことはわからない。混入する度合いは実験ごとに異なり、キットの改良によって次第に下がっていくだろう。定性的に、こういう転写産物があるということを知ることしかできない。それでも、オルガネラゲノムの配列が知られていない生物種の場合はとても役に立つだろう。またモデル生物やヒトでもオルガネラゲノムに生じている変異を検出するには十分使えるのではないか(定量性はないけれど)。

上に書いた Smith DR 博士らによって、オルガネラ転写産物を分析するためのソフトウェア Chloroseq が開発された。   http://github.com/BenoitCastandet/chloroseq   論文としても発表されている。   https://www.ncbi.nlm.nih.gov/pubmed?linkname=pubmed_pubmed&from_uid=27402360    ChloroSeq, an Optimized Chloroplast RNA-Seq Bioinformatic Pipeline, Reveals Remodeling of the Organellar Transcriptome Under Heat Stress.   Castandet B, Hotto AM, Strickler SR, Stern DB.   G3 (Bethesda). 2016 Sep 8;6(9):2817-27. doi: 10.1534/g3.116.030783.   PMID: 27402360

RNA-seq ではサンプルにゲノム DNA が残存するとか、サンプル自体に微生物が混ざっていてそれらに由来する配列が出てくるとか、主目的でない要因で出てくる配列が混ざっていることがあるらしい。画像や音声のデータでは、混ざったデータを分離する解析法が開発されている。

何か巧妙な計算によって、オルガネラゲノムから転写される mRNA に関する定量的な情報を得られるようにできればオルガネラ研究、また核とオルガネラの相互作用の研究にとても役に立つと思うが、無理なようにも思える。 何が問題か: オルガネラ由来の mRNA がサンプルに混入する度合いは実験ごとに異なると推測される。また mRNA の種類が異なると変化してしまうだろう。 オリゴ dT をプライマーにしている場合、ポリA がない RNA では A が連続している部分からのみ DNA 合成がおきて、それ以外の場所からはおきにくいという偏りが生じる。 なにか「基準になる RNA・値」を決めておいて相対的な比率を見るとしても、どんな基準があるか。全く望みがないわけでもない(信用出来ない値でも、ないよりはましと考えれば: p-value のような信頼度を表す値を付けられればよい)かもしれない。

ある経済学者のいうことによると、こういうことを考えるには「単純化してモデルを扱い易くするための仮定」を見つける才能が必要らしい。それは「そんなことは実際にはあり得ないじゃないか」といわれるようなものでも全然問題ない。例えば「オルガネラ由来の mRNA がサンプルに混入する度合いは、どの種類のオルガネラ mRNA でも同じだ」とかが考えられる。

また、現在では様々な non-coding RNA の機能が解明されたことによりポリA 鎖がついていない RNA の解析の重要性が増している。そのためポリA 鎖をもつ RNA 以外も標的にしたデータも増えている。

例: E-GEOD-23439 - Genome-wide double-stranded RNA sequencing, E-GEOD-49325 - Characterization of stress-responsive lncRNAs というデータがあった。

High-throughput m6A-seq reveals RNA m6A methylation patterns in the chloroplast and mitochondria transcriptomes of Arabidopsis thaliana.   Wang Z, Tang K, Zhang D, Wan Y, Wen Y, Lu Q, Wang L.   PLoS One. 2017 Nov 13;12(11):e0185612. doi: 10.1371/journal.pone.0185612. eCollection 2017.   PMID: 29131848 という論文では、GSE72706 という公開データを分析している。

ArrayExpress では、「Type」の部分をクリックするとタイプで並び替えられて、「RNA-seq of non coding RNA」のタイプのデータを集めることができる。miRNA を分析しているデータが多い。


SRR929000 のサンプルは、ZFP3 overexpressing and Col-0 wild type seeds were germinated on media supplemented with 5uM estradiol or 5uM estradiol and 2.5uM ABA.   エストラジオール誘導性プロモーターを使って ZEP3 発現を ON, OFF している。「ZEP3 が強く発現している状態に、ABA 処理したもの」と「ABA 処理のみしたもの」を比較することになる。  エストラジオールだけを与えて育成したものも作っている。それは成長が早いので 36時間後にサンプリングしている。 「ZEP3 が強く発現している状態に、ABA 処理したもの」と「ABA 処理のみしたもの」はABAのせいで生育が遅いので 84 時間後にサンプリングしている。

ArrayExpress には本体データはなく、リンクされている GEO http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48661 に、SRA ファイルへのリンクがある。 これをダウンロードするには、SRA Toolkit という専用のソフトウェアを使う。SRR929000, SRR929001, SRR929002, SRR929003 という番号が、4つのファイルにそれぞれついている。これらの番号をメモしておく。929000 は Col-0/control 929001 は Col-0/ABA 929002 は ZFP3ox/control 929003 は ZFP3ox/ABA

sra ファイルとはどんなものか。http://www.ncbi.nlm.nih.gov/books/NBK47540/   が説明である。Sequence Read Archive を略して SRA という。 http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software   SRA Toolkit などのソフトウェアが公開されている。これらのものを開発して頂いているおかげで、データを誰でも操作できるようになっている。

Run Browser というページ http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=run_browser Web browser でデータを見ることもできる。このページで、今から分析したい SRR929000 を入力すると、内容のページに飛ぶ。

前処理、分析のためのソフトウェア

SRA Toolkit

SRA Toolkit をインストールしてみる。

VirtualBox というソフトウェアに、BioLinux という Linux のディストリビューションを入れて使うこともできる。これは RNA-seq などを分析するためのプログラムが最初から入っている。しかしそれらのソフトウェアを管理することが難しい。

BIOCONDA https://bioconda.github.io/ という便利なソフトウェアセンターが公開されている。 ソフトウェアを一つずつばらばらに入れていると管理しにくい。スマートフォンでアプリを入れるのと同じように、どんなソフトウェアもソフトウェアセンターから入れる方が管理しやすい。 BIOCONDA は、Miniconda, Anaconda という管理システムを使って、生物学関連のソフトウェアを集めて簡単にインストールできるようになっている。

Ubuntu 16.04 をインストールした PC を用意した。それに BIOCONDA から生物学関連のソフトウェアをインストールする。こちらの方が、最新のソフトウェアに簡単に確実に更新することができる。 しかし、tophat というソフトウェアを使おうとしたら、BIOCONDA に入っていたものではうまく動かず、tophat のホームページから取ってきたもので動いた(python2 で動かすように修正が必要だった)。その原因は、Python3 用の conda を入れたせいだった。

Miniconda には二種類あって、「Python2 用の Miniconda」をインストールすれば修正しなくてもよくなる。このことについては、この文章よりもっと詳細でわかりやすい解説が公開されている。   http://imamachi-n.hatenablog.com/entry/2017/01/14/212719 「biocondaを利用してNGS関連のソフトウェアを一括でインストールする」 Imamachi-n 氏による解説

また、Imamachi-n 氏の解説にあるように、Python のバージョンが異なる仮想環境を作ることもできる。私は Imamachi-n 氏の解説を見て、以下のようにした。

 $ conda create --name py27 python=2.7
 py27に入る     $ source activate py27
 py27から抜ける $ source deactivate py27

[py27] がプロンプトに出た状態で「conda install tophat」とする。py27 を activate した状態では、Python2.7 を使う環境で tophat が動くようになり、エラーは出なくなった。

ソフトウェアをアップデートするときは、「conda update sra-tools」のようにコマンドを打って行う。

SRA Toolkit は、コマンドラインから使うプログラムの詰め合わせである。これも BIOCONDA に入っていた。

http://www.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=toolkit_doc&f=std#s-3   SRA Toolkit Installation and Configuration guide にあるように、テストしてみる。

 fastq-dump -X 5 -Z SRR390728

と打つと、データの5行分が表示された。今から分析したい SRR929000 でもうまくいった。このように一部を見るだけなら待ち時間は少ない。この方法で出てくるデータは、fastq というフォーマットで記述されている。細かい型式はシーケンサーが異なる会社の機種だと変わる。 SRA Toolkit には、機種それぞれに適合させたソフトウェアも用意されている。 SOLID のデータの場合、abi-dump というツールが用意されている。

 abi-dump -v SRR929000

と打つと、SRR929000_F3.csfasta というファイルと、SRR929000_F3_QV.qual という2つのファイルがダウンロードされた。.qual の方は quality data を含んでいる。 SOLID のデータを分析する場合、この 2 つのデータを必要とする場合がある。

fastq について

一つの read は、4行で表現される。一行目は、@ から始まり、説明文になっている。次の行は塩基配列である。しかし、SOLID の場合、ATGC ではなく T021100030..... というような数字の列になっている。これは塩基の読み方のせいである。 SOLID に対応しているプログラム (bowtie) もあるし、していないものもある。 3行目は + で始まり、1行目と同じことが書いてある。 4行目は、データ(塩基)の信頼性を表す情報が入っている。目で見ただけでは何のことか分からない。

いろいろなファイル

(multi-) fasta 形式ファイルは、大量の塩基配列とそれらの簡単な注釈(一行ずつ)を含むデータで、様々な分析に供せられる。

fastq 形式ファイルは、大量の塩基配列とそれらの信頼性とそれらの簡単な注釈(一行ずつ)を含むデータで、様々な分析に供せられる。

SRAファイルのアクセッション番号が分かれば、SRA Toolkit の fastq-dump を用いて、直接 fastq ファイルを取り込める。

SOLID のデータは、abi-dump で取り込むこともできる。csfasta ファイルと、qual ファイルの二つが取り込まれる。csfasta は ABI のマシンの配列データで、qual ファイルは品質を示している。

fai ファイルは、fasta ファイルから生成したインデックスを含む。染色体の名前、染色体の長さ、そのデータのスタート位置(ファイルの先頭から数える)、一行の文字数、一行のバイト数(文字+改行コードで数える) が、データごとに書かれている。

SAM 型式ファイルは、bowtie, bowtie2, tophat, hisat2 のようなマッピングソフトウェアの出力である。それをさらに BAM 形式に変換して、カウント情報の取得などを行うことが多い。SAMtools http://samtools.sourceforge.net/ というソフトウェアで変換できる。 SAM ファイルはテキストデータになっていて、人間が一応読むことができる。BAM ファイルは SAM ファイルと同じ情報を持っているが、ソフトウェアで処理しやすいバイナリファイルになっていて、人間が見てもわからない。SAMtools は、BAM ファイルが持っている情報の必要な部分だけを出力するためにも用いられる。その出力をエクセルや他のプログラムで処理することもできる。

BED 形式ファイルは、ゲノム上の特定の領域を示す情報が書かれているファイルで、転写因子の結合部位などを示すのに使われる。BAM ファイルと対になっていることが多い。ChIP 実験の結果のファイルには、BAM ファイルと BED ファイルが同梱されていることがある。配列データを tophat で処理すると、accepted_hits.bam, unmapped.bam ファイルと同時に deletions.bed, insertions,bed, junctions.bed ファイルが生成する。

NGS Surfer's Wiki https://cell-innovation.nig.ac.jp/wiki/tiki-index.php   で、詳しくわかりやすく解説されている。

GSM861508_PM1_m1_btb_chrom.bed というファイルでは、8601636 行のデータが納められていた。目で眺めても全く意味をつかむことはできない。BED 形式ファイルは、他のプログラムの入力ファイルとして使われると NGS Surfer's Wiki に書かれている。 他のプログラムによって、人間が見てわかるような形式で表示する。

GFF/GTF形式ファイル は、遺伝子アノテーションの情報を含む。http://ccb.jhu.edu/software/tophat/index.shtml から、「Index and annotation downloads」へ進んで、それぞれの生物用のファイルを入手する。これは Cuffdiff というプログラムの入力の一つとして使われる。

refFlat形式ファイルは、遺伝子アノテーションの情報を含む。

fastq-dump のオプションについて

fastq-dump --help で見ると、たくさんのオプションがある。 例:

 fastq-dump -X 5 -Z SRR390728

最後にアクセッション番号を書いて、取り出したいデータを指定する。

とりあえず、-X 10000 として小さめの fastq ファイルを作ってみる。3秒くらいで作れた。サイズは 2.4Mbyte になった。

データ全体(1660万のリードを含む)を処理すると、約2時間で終了して 4GB の fastq ファイルになった。

 >fastq-dump SRR929000
 Read 16603126 spots for SRR929000
 Written 16603126 spots for SRR929000

マッピング

様々なプログラムで、RNAseq のデータをゲノム配列にマッピングできる。 Bowtie というプログラムが代表的なものである。 しかし、スプライシングアイソフォームのことを考えに入れる場合は、それだけでは不十分であると解説されている。RNAseq ではスプライシングの違いを検出するための情報が得られるので、それを生かさないのはもったいない。そこでそういうスプライシングに対応したプログラムが開発されている。

マッピングプログラムには、二種類の入力ファイルを与える。一つは RNAseq データで、もう一つはマッピングされる側のゲノム配列である。

マッピングされる側のデータは、どういうものか。ゲノム全体配列が入っているデータを元にして、それをインデックスファイルにしたものを使う。

インデックスファイルは、マッピングプログラムごとに別のものを用意する必要がある。 bowtie と bowtie2 は、名前は同じようだが全く異なるもので、インデックスファイルには互換性が全然ない。 それぞれのためにインデックスファイルを作らないといけない。

Bowtie について

Bowtie は短いリードを含むデータ(初期のマシンからのデータに多い)に適していて、SOLID のデータも処理できるので、これをまず使ってみる。

Bowtie のホームページ   http://bowtie-bio.sourceforge.net/index.shtml

上に書いたように BIOCONDA を使えるようにしてあれば、

 conda install bowtie

BIOCONDA から導入した場合は、サンプルファイルが得られないという問題点がある。

Bowtie を実行するには、入力ファイルを用意する必要がある。一式を手動で導入した場合、サンプルファイルが同封されているので説明にあるように以下のコマンドを実行してテストする。

 bowtie e_coli reads/e_coli_1000.fq | less

e_coli で、大腸菌ゲノム配列をインデックスファイルにしたものを使うことを指定する。indexes ディレクトリに、e_coli で名前が始まる ebwt ファイルが6つある。これは大腸菌ゲノム全体が、インデックスファイルにされている。

reads/e_coli_1000.fq で、reads ディレクトリにある fastq ファイルを指定する。1000 個のリードの配列が ATGC..... で記載されている。

そうすると、すぐに Bowtie 型式の出力が出てくる。最初の行は、

 r0  (Tab)  -  (Tab) gi|110640213|ref|NC_008253.1|  (Tab)  3658049  (Tab) ATGCTGGAATGGCGATAGTTGGGT
 GGGTATCGTTC   (Tab)   45567778999:9;;<===>?@@@@AAAABCCCDE   (Tab)   0  (Tab)   32:T>G,34:G>A

となっていた。

タブ区切りのテキストデータになっている。最初から 「そのリードの注釈(この場合注釈がないので0番目のリードで r0)、マッチした strand (+ or -)、インデックスファイルの元になったデータの名前、マッチした場所、塩基配列、品質を表す記号、信頼できない塩基の数?、ミスマッチがある場合その報告」 となっている。一番最初のリードは、fastq ファイルを見ると GAACGATACCCACCCAACTATCGCCATTCCAGCAT という配列である。それの相補鎖が、結果に書かれている。その前に  3658049 と書かれているのは、 3658049 番目の塩基の部分から、ゲノムと最初のリードがマッチしているということを表している。

オプションに --sam とつけると、出力が Sequence Alignment/Map (SAM) format に変わる。 このフォーマット、それを他のソフトウェアで使いやすいバイナリファイルに変換した BAM がよく使われる。

SAM フォーマットの出力の最初の部分:

 @HD     VN:1.0  SO:unsorted
 @SQ     SN:gi|110640213|ref|NC_008253.1|        LN:4938920
 @PG     ID:Bowtie       VN:1.0.1        CL:"bowtie e_coli reads/e_coli_1000.fq --sam"
 r0      16      gi|110640213|ref|NC_008253.1|   3658050 255     35M     *
 0       0       ATGCTGGAATGGCGATAGTTGGGTGGGTATCGTTC     45567778999:9;;<===>?@@@
 @AAAABCCCDE     XA:i:0  MD:Z:0G1T32     NM:i:2

@ から始まるのはヘッダーの部分で、今回使っている大腸菌ゲノムのデータのアクセッション番号が出ている。Bowtie で処理されたことが記載されている。 その後にマッピングの結果が出ている。

つぎに、実際のデータでテストする。

bowtie のためのインデックスファイルは、bowtie-build というプログラムでできる。試しにやってみた。オルガネラゲノムだと小さいのですぐに処理できた。モデル生物なら、すでにゲノム全体がインデックスファイルにされたものを Bowtie のページから入手できる。

SRR929000 のデータは SOLID のデータなので、普通のインデックスされた参照配列ではだめで、カラースペース用のインデックスファイルでないといけない。 Fasta ファイルがあれば、それを元に bowtie-build  ゲノムファイル名  インデックスのベースになる名前  で、インデックスファイルを作ることができる。カラースペース用ではオプションに -C をつける。小さいゲノムデータとしてはオルガネラゲノムが使える。それなら一瞬で処理が完了する。できた ebwt ファイルを適当なディレクトリに移動する。

https://github.com/brentp/bowfast/tree/master/aligner-compare   に、bowtie のオプションの付け方について解説がある。それを参考にしてみる。

 #!/bin/bash
 SRR=SRR929000
 bowtie   -C   /インデックスファイルのある場所/インデックスファイルのベースになる名前    -Q  /ファイルのある場所/${SRR}_F3_QV.qual    -f /ファイルのある場所/${SRR}_F3.cfasta    --chunkmbs 1025   --best   --strata  -v 1  -m 1  -y  -n 1  -p 4  --maxbts 25000   /出力ファイルの場所/${SRR}_名前

SOLID のデータの場合 -C を付ける。それ以外の機種のデータの場合は bowtie 以外のソフトウェアを使うほうがよいらしい。bowtie は短いリードを処理するのに向いていると書かれていた。

インデックスファイルがオルガネラゲノム由来で小さいと、すぐに処理が完了した。ゲノム全体を処理するならそれだけ時間がかかる。 実際にやってみるとオルガネラゲノムにマッピングされるリードの割合はとても少ない。またマッピングされる場所に偏りがある。 これはとても残念なことだが、仕方がない。

Tophat

Imamachi-n 氏の解説にあるように、conda を用いて Python のバージョンが 2.7 の仮想環境を作る。そこに tophat をインストールする。元の conda が Python3 の場合そうしないとエラーが出た。

 $ conda create --name py27 python=2.7
 py27に入る     $ source activate py27
 py27から抜ける $ source deactivate py27

[py27] がプロンプトに出た状態で「conda install tophat」とする。py27 を activate すると、Python2.7 を使う環境で tophat が動くようになり、エラーは出なくなった。

$ source activate py27 として、[py27] がプロンプトに出た状態にする。そしてコマンドを打つ。

 tophat -o 出力するディレクトリー  インデックスファイルの名前   入力するファイルの名前

この場合のインデックスファイルは、bowtie2-build で作成する。

指定したディレクトリーに、複数のファイルができる。accepted_hits.bam が、一番大切な結果のファイルになる。 bam ファイルから人間が読める形で情報を取り出すには、samtools というコマンドを使う。

SRR064979 というデータを tophat で処理してみた。このデータは Illumina genome analyzer のデータで、E-GEOD-23439 - Genome-wide double-stranded RNA sequencing というタイトルがついている。 リボソーム RNA を除去した後に、RNA の一本鎖の部分のみを分解する RNase One という RNase を用いて処理している。 このデータでは、オルガネラゲノムにマッチするリードはかなり多かった。10 分くらいでマッピングは完了した。

HISAT2

この場合も、インデックスファイルを用意する。hisat2-build というコマンドを使う。オルガネラゲノムの fasta ファイルを用いて試してみる。

 hisat2-build -f /ファイルの場所/ファイル名.fasta /出力するディレクトリー/インデックス名

Manual に、次のように書かれている。

 hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <hit>]

入力ファイルがペアになっている場合、-1 と -2 で指定する。そうでない場合 -U で指定する。 または、--sra-acc の後ろに SRA の accession number を書いてもよい。-S で、出力となる SAM ファイルの場所と名前を指定する。

そうすると、一晩経っても計算が終了しなかった。Ctrl-C で止めると結果として巨大な SAM ファイルが作られていた。どこが悪いのかというと、単にデータが巨大なせいらしい。

このデータを fastq-dump で取り込もうとした。丸 2 日かかっても終わらない。大きすぎて分析に向いていない。こういうデータもあるので注意しないといけない。

samtools の使い方

マッピングソフトウェアの出力である sam, bam ファイルから人間が読める形で情報を取り出すには、samtools を使う。

cufflinks

これについては、AtCAST3:Home   http://atpbsmd.yokohama-cu.ac.jp/cgi/atcast/home.cgi   植物の遺伝子解析のための、横浜市立大学のすばらしいページの「Data format for RNA-seq data」のセクションで解説されている。

マイクロアレイデータでは、「遺伝子1のID、条件 A での発現数値、条件 B での発現数値、」で一行が構成され、それが遺伝子の数だけ並んでいる。 RNA-seq データも、同じようなフォーマットにできると便利である。エクセルにも入るようになる。AtCAST3 でも、そういうフォーマットにすると分析できる。 「Data format for RNA-seq data」のセクションに、「Tutorial: Calculation of Gene counts using cufflinks」がある。これに従えばよい。

Bowtie や SAMtools の他に、tophat2、 cufflinks をインストールする。

Arabidopsis ENSEMBL TAIR10 genome and annotation files を http://ccb.jhu.edu/software/tophat/index.shtml  から入手する。圧縮ファイルなので元に戻して、適当に配置する。

SRA、または fastq ファイルを用意する。わかりやすいように名前を変える。

例として掲載されているコマンドを参考にして処理する。

Fastq ファイルを tophat で処理する。これで accepted_hits.bam ファイルができる。

BAM ファイルを cuffdiff コマンドで処理する。生成したファイルは、「遺伝子1のID、条件 A でのカウント数値、条件 B でのカウント数値、」 になっている。こうできれば、マイクロアレイデータに近くなる。

BAM ファイルとはどんなものか。SAM 型式ファイルは、マッピングソフトウェアの出力である。それをさらに BAM 形式に変換して、カウント情報の取得などを行うことが多い。BAM ファイルと SAM ファイルは同じ情報を含んでいるが、BAM ファイルはソフトウェアで処理しやすい形になっている。

accepted_hits.bam ファイルを cuffdiff コマンドで処理する。生成したファイルは、「遺伝子1のID、条件 A でのカウント数値、条件 B でのカウント数値、」 になっている。こうできれば、マイクロアレイデータに近くなる。 cufflinks の使い方については、AtCAST3:Home   http://atpbsmd.yokohama-cu.ac.jp/cgi/atcast/home.cgi   植物の遺伝子解析のための、横浜市立大学のすばらしいページの「Data format for RNA-seq data」のセクションで解説されている。 変換の流れは、http://cole-trapnell-lab.github.io/cufflinks/getting_started/   cufflinks のホームページで解説されている。

すでに accepted_hits.bam ファイルが用意されている場合、「Data format for RNA-seq data」の下の方に例としてこう書かれている。

 cuffdiff ./Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Annotation/Genes/genes.gtf tophat_wt_rep1/accepted_hits.bam  tophat_wt_rep2/accepted_hits.bam  tophat_wt_rep3/accepted_hits.bam  tophat_pifq_rep1/accepted_hits.bam  tophat_pifq_rep2/accepted_hits.bam  tophat_pifq_rep3/accepted_hits.bam &&

cuffdiff コマンドは、一つの引数として gtf ファイルを取る。このファイルは遺伝子アノテーションの情報を含む。それ以外に、発現量の情報を含む複数の bam ファイルを引数に取れる。上の例では二種類のサンプルについて、それぞれ 3 回の実験結果の bam ファイルがあるので 6 つの bam ファイルを一度に扱っている。これらの bam ファイルはマッピングソフトウェアの出力だから、リードのカウント数(発現量の情報)をそれぞれの遺伝子に紐付けしてくれることになる。

シロイヌナズナの場合、gtf ファイルは、Arabidopsis ENSEMBL TAIR10 genome and annotation files を http://ccb.jhu.edu/software/tophat/index.shtml  から入手する。圧縮ファイルなので元に戻す。genes.gtf ファイル以外にもたくさんの種類のファイルが入っている。 普通に中身をコピーすれば解説にあるように、./Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Annotation/Genes/genes.gtf に配置される。

マイクロアレイデータでは、RMA 法のような前処理(標準化とノイズ除去)法がある。RNA-seq ではどうなのか、これから調べる。

(続く)

vim: set ts=8 sts=2 sw=2 et ft=a111_modified_flexwiki textwidth=0 lsp=12: