RNA-Seqの数理―生成モデルによる発現量推定：アーカイブ

この「RNA-Seqの数理」シリーズでは，次世代シーケンサを用いたRNA-Seqにおける発現量推定の数理モデルを理解することを目的とする．

副題にある「生成モデル」とは，観測データの生成過程を確率的にモデル化し，データが与えられたときの事後確率を用いて分類したいクラスや予測したい値を推定する方法のことを指す．今回のRNA-Seqにおける生成モデルでは，次世代シーケンサで読み取られた配列が，どのようにして細胞内のトランスクリプトームから実験によって読み取られて観測されたかという一連の流れを，生成モデルとして表現する．そして確率的なもっともらしさやパラメータの推定をおこない，トランスクリプトームの発現量を求める．

このシリーズについて

次世代シーケンサによるRNA-Seqの発現量推定といっても，実験対象は大腸菌レベルから人に至るまでゲノムサイズや遺伝子数は多種多様であり，実験機器も各メーカーごとに独自に改良が加えられてきた．しかしながら，RNA-Seqはトランスクリプトームの配列を超並列に大量に解読するという点で，基本的な原理は生物種やプラットフォームに依存しない．まずはこの共通部分に関して単純なモデルを作成し，そこから個別の生物の特性や実験方法に合わせた複雑なモデルを構築していく．

今回のシリーズは前提条件として，発現量推定の対象とする種をゲノム既知なモデル生物に限定する．つまり，ヒトやマウス，シロイヌナズナ，線虫など，すでにゲノム配列が解読されていて，なおかつ遺伝子の配列や機能がよく研究されている種を対象にするということだ．この仮定をおかないと，そもそもデータ処理の時点でショートリードのマッピングができなかったり，モデルが立てられないということになりかねない．とはいうものの，仮定を立てたからといってマッピングの曖昧さは解消されないし，アイソフォームの存在や選択的スプライシングなど複雑な転写メカニズムはいくらでも考えられるので，それらは追々モデルに組み込んでいく．

注意

このシリーズの記事で書かれていることはすべて個人的な勉強記録であって，数理的な解釈や導出の厳密性を保証するものではない．勘違いや間違いなどが多分に含まれている可能性があるので，参考にされる際には十分に注意していただきたい．もし間違いや誤字脱字等見つけられた方は，コメント欄などで指摘していただければ幸い．

このシリーズの目標

このシリーズの目標（そして個人的な勉強の目標）としては，

  • CufflinksやRSEMなどのスタンダードとなっている発現量推定の数理的背景を理解する
  • 発現差異解析(DE)などの複数サンプルにまたがる発現量比較の手法の理解につなげる

ことを考えている．基本的には論文に記載されている内容をベースに進めていくので，論文のmethodの理解が最終的な目標となるだろう．

記事まとめ

参考



ディレクトリの書き込み権限が無いとファイルの消去ができない

久しぶりに初歩的なところでハマったのでメモ．多人数が使用するサーバ上で，何故かrmすることのできないファイルがあった．パーミッションも所有者もグループも特に変なところは見当たらず，ファイルが壊れているわけでもGUIで他のプログラムが使用中というわけでもない．なのに，rmしようとすると以下のようなエラーが出てしまう（例示しているコマンドラインは手元で再現した状況）．

1
2
3

$ rm hehehe
remove hehehe? y
rm: hehehe: Permission denied

sudo権限を持っていないのでなんともしがたいなーと思いつつ，適当にディレクトリ名を”fuga_broken”に変更して放置していたのだけれども，たまたま機会があって人に相談したら「ディレクトリの書き込み権限が無いんじゃない？」と言われて，そこでようやく気がついた．

1
2
3

$ ls -l
合計 0
dr-xr-xr-x 3 yag_ays staff 102  7 13 15:27 fuga/

というわけで，ディレクトリの書き込み権限がなかったせいで，そのディレクトリ内のファイルの消去ができなかっただけだったというオチ．なぜ書き込み権限が無くなっていたのかはいまだに不明だ．対象となるディレクトリをchmodで書き込み権限を付加したら，あっさりファイルを削除することができた．

1
2
3
4

$ chmod +w fuga
$ rm fuga/hehehe
remove fuga/hehehe? y
$

なんともまあ呆気なく解決してしまった．普段からディレクトリのパーミッションはまったく変えないので，問題解決の際にディレクトリにまで考えが及ばなく，久々に自分の経験の浅さを実感した．



Tophat/Cufflinksなどの開発者であるCole TrapnellによるRNA-Seqの遺伝子発現解析に関する講義動画

非常に良い講義資料およびその講義ビデオを見つけたので紹介する．これはハーバードの大学院およびハーバード公衆衛生大学院のComputational Biologyの講義で，マイクロアレイ/次世代シーケンサの発現解析に始まり，ChIP-Seqやモチーフ検索，GWAやSNP解析に至るまで，Computational Biologyにおける様々な研究対象に関する講義が行われている．

STAT115 Introduction to Computational Biology

This course is an introduction to state-of-the-arts concepts, methodologies and practices in computational biology. After taking this course, you should be able to:

  • Use appropriate commands or write your own scripts to manipulate your own genomic data

  • Exploit appropriate pipeline or software to analyze your own high-throughput gene expression microarray and/or RNA-Seq data and high-throughput sequencing transcriptional regulation data

  • Apply statistical or mathematical models to your own population/human genetics data

  • Criticize and/or extend the computational analysis or statistical models in the literature

  • Use R, Python, and/or Shell scripts and basic concepts in genomics and genetics to perform your own bioinformatics and computational biology research

http://www.stat115.org/syllabus.html

このなかで，Cole TrapnellがRNA-Seqの遺伝子発現解析に関してレクチャーしている．Cole Trapnellといえば，TopHatやCufflinksの開発者の中心人物のひとりだ．これらのツールに関する論文の第一著者でもある．その彼が，RNA-Seqの定量化や発現量推定などの一連の流れを，自身が開発に参加した”Tuxedo Tools”（Bowtie, Tophat, Cufflinks, CummeRbundの総称）を例に紹介しているのがこの動画だ．主にマッピング後の発現量推定に関して，基本的なFPKMの求め方であったり，Cufflinksのisoformの取り扱いであったりを，少し数式も交えながら詳しく解説している．

講義資料：http://www.stat115.org/lectures/stat115_rnaseq.pdf

lecture9-I from stat115 on Vimeo.

lecture9-II from stat115 on Vimeo.



TimeTreeで生物種間の分岐年代を調べる

2種類の生物がどれくらい系統的に近いか遠いかといった情報は，生物の分類をWikipediaで確認したり系統樹を眺めたりと色々と確認する方法があるが，ではどれくらい前に2種類の生物が分岐したかという情報は意外と確認するのが難しい．そういった分岐年代を簡単に確認，比較することができるのがTimeTreeというサイトだ．

TimeTree :: The Timescale of Life

生物種を入力する

TimeTreeウェブサイトの左上にある検索ボックスに，分岐年代を辿りたい生物種を入力する．通常はCanis lupusとFelis catusのような学名を入力するのだけれども，一般的な名称であるdogやcatも自動で変換してくれるようだ．もし学名がハッキリ分かっているのなら，そちらを入力したほうが確実だろう．

結果

今回のような検索クエリを入力した結果，イヌ(Canis lupus familiaris)とネコ(Felis catus)の分岐年代は55.1 Million Years Agoと推定された．ただしこの値というのは，様々な分岐年代の研究に関してイヌとネコが当てはまるものを複数集めてきて，最終的に平均を取った値となっている．個別の研究の推定結果というのは，図の左側の時系列の箇所に黒の点で表示されているように，実際には42Myaから67Myaまで幅があるようだ．このように複数のソースから分岐年代を確認できるのもTimeTreeの大きな特徴となっている．

上図の黒の点で示した個別の研究というのは，下のリストからそれぞれ確認することができる．Taxon Aがイヌ，Taxon Bがネコの検索クエリにそれぞれ対応しており，研究によってはイヌ亜目Caniformiaやネコ亜目Feliformiaのような目レベルで確認されているものや，クマ科Ursidaeやアザラシ科Phocidaeといった少し系統的に離れた分類群との比較も列挙されている（クマやアザラシはどちらもネコ目）．

このように，ウェブページに表示されているリストから詳細なページに飛ぶことができるほか，「Molecular Time Estimates」にある「Download as CSV」をクリックすると，根拠となる論文や出版年のリストをコンマ区切りテキストとしてダウンロードすることができる．

参考



ゲノムアセンブリにおいて最適なk-merを推定するKmerGenieを試してみた

KmerGenieという，シーケンスデータのみを用いてアセンブラに依存しない形でアセンブルに指定する最適なk-merの値を求めることができるソフトウェアが最近出てきたので，ちょっと使ってみた．

KmerGenie

そういえば以前VelvetKという，これも同じようにk-merパラメータをデータセットから推定するツールについて書いたが，これは名前の通りVelvetに特化したソフトウェアだった．2007年くらいに公開され使われ始めたVelvetは，NGSの普及と相まってかなり長い間ゲノムアセンブラとして広く使われてきた印象がある．その後ABySSやMIRA，ALLPATHS-LG，Rayなどの様々なゲノムアセンブラが出てきたが，それらに対してVelvetKのようなk-mer最適化のようなツールは無かったようで，そのあたりアセンブラに非依存というのがKmerGenieの売りのようだ．

KmerGenieの仕組み

KmerGenieがどういう仕組みで最適なk-merを推定しているかは，以下のスライドやarXiv.orgで公開されている論文に詳しく書かれている．

これによると，単純に言ってしまえば最適なk-merの出現頻度というものはきれいな正規分布に従うと仮定した上で，k-merの値を色々変えて出現頻度のヒストグラムを作成し，生成モデルを立ててフィッティングをして最適なk-merの値を推定しているようだ．実際には，正規分布を仮定したゲノム配列の分布とパレート分布を仮定したエラーの配列の分布が混ざった混合分布を考えて，その混合の重みなども考慮したパラメータ最適化をしているらしい．

KmerGenieのインストール

まずはKmerGenieをインストールする．ダウンロードしたファイルを展開すると既にディレクトリにkmergenieという実行ファイルがあるが，makeをしないと利用することができないので注意．

1
2
3

$ wget http://kmergenie.bx.psu.edu/kmergenie-1.5397.tar.gz
$ cd kmergenie-1.5397
$ make

KmerGenieの実行

それでは実際にfastqファイルに対してKmerGenieを実行してみる．

1. ゲノムアセンブリのチュートリアルで使用されたE. coliのゲノムシーケンス

まずは2013年の6月に行われたMSU NGS Summer Course 2013のチュートリアルで使用されたサンプルデータを使ってみる．これはE. coliのゲノムシーケンスで，元はChitsaz et al.が出したデータの一部のようだ．中身は約4.7Mの70bpのシングルリード．

Assembling E. coli sequences with Velvet — ANGUS 2.0 documentation

1
2

$ curl -O https://s3.amazonaws.com/public.ged.msu.edu/ecoli_ref-5m-trim.fastq.gz
$ gunzip ecoli_ref-5m-trim.fastq.gz

それではecoli_ref-5m-trim.fastqに対してKmerGenieを実行する．実行後に作られる各k-merでのカウントデータや分布をプロットしたpdfはカレントディレクトリにそのまま出力されるので，もし必要ならディレクトリを別途作ってそこで実行したほうが良いだろう．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

$ kmergenie ecoli_ref-5m-trim.fastq
running histogram estimation
Linear estimation: ~130 M distinct 41-mers are in the reads
K-mer sampling: 1/100
| processing                                                                                         |
[going to estimate histograms for values of k: 61 51 41 31 21 
-----------------------------------------------------------------------------------------------------------------------------Total time Wallclock  77.7066 s
fitting model to histograms to estimate best k
fitting histogram for k = 21
fitting histogram for k = 31
fitting histogram for k = 41
fitting histogram for k = 51
fitting histogram for k = 61
estimation of the best k so far: 51
refining estimation around [45; 57], with a step of 2
running histogram estimation
Linear estimation: ~139 M distinct 39-mers are in the reads
K-mer sampling: 1/100
| processing                                                                                         |
[going to estimate histograms for values of k: 57 55 53 51 49 47 45 
-----------------------------------------------------------------------------------------------------------------------------Total time Wallclock  56.4315 s
fitting model to histograms to estimate best k
fitting histogram for k = 21
fitting histogram for k = 31
fitting histogram for k = 41
fitting histogram for k = 45
fitting histogram for k = 47
fitting histogram for k = 49
fitting histogram for k = 51
fitting histogram for k = 53
fitting histogram for k = 55
fitting histogram for k = 57
fitting histogram for k = 61
table of predicted num. of genomic k-mers: histograms.dat
best k: 55

今回の場合，KmerGenieによる最適なk-merの値は55となった．出力結果で示されるk-mer頻度のプロットを並べて見ると，確かにkが小さいときは少しいびつになっており，k=55あたりで一番正規分布っぽくなっている…??ことがわかる．

ただ，このサンプルーデータを解析しているチュートリアルではk-merを31〜35にしてvelvetでアセンブルしていることから，今回の55という推定値は少し大きめのような気がする．

2. Assemblathon2で使用されたセキセイインコのゲノムシーケンス

次はもう少し大きいデータで試してみよう．Assemblathon2というアセンブリツールの性能を競うコンペで使用されたセキセイインコのゲノムシーケンスの中からDuke Illumina GAIIx runsというデータを使ってみる．中身は76bpのペアエンドで合計で126M readsある．

Index of /Data/hcbxz0i7kg/Parrot/illumina_duke_runs

1
2

$ lftp -e "mirror" http://bioshare.bioinformatics.ucdavis.edu/Data/hcbxz0i7kg/Parrot/illumina_duke_runs/ 
$ cat *.fastq > budgie.fastq

Assemblathonのサイトによるとセキセイインコのゲノムは二倍体らしいので，今回はKmerGenieの「–diploid」のパラメータを付けて実行する．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

$ kmergenie --diploid budgie.fastq
running histogram estimation
Linear estimation: ~3721 M distinct 46-mers are in the reads
K-mer sampling: 1/1000
| processing                                                                                         |
[going to estimate histograms for values of k: 71 61 51 41 31 21
---------------------------------------------------------------------------------------------------------------------------Total time Wallclock  3035.49 s
fitting model to histograms to estimate best k
fitting histogram for k = 21
fitting histogram for k = 31
fitting histogram for k = 41
fitting histogram for k = 51
fitting histogram for k = 61
fitting histogram for k = 71
estimation of the best k so far: 21
refining estimation around [15; 27], with a step of 2
running histogram estimation
Linear estimation: ~6335 M distinct 24-mers are in the reads
K-mer sampling: 1/1000
| processing                                                                                         |
[going to estimate histograms for values of k: 27 25 23 21 19 17 15
---------------------------------------------------------------------------------------------------------------------------Total time Wallclock  2305.61 s
fitting model to histograms to estimate best k
fitting histogram for k = 15
fitting histogram for k = 17
fitting histogram for k = 19
fitting histogram for k = 21
fitting histogram for k = 23
fitting histogram for k = 25
fitting histogram for k = 27
fitting histogram for k = 31
fitting histogram for k = 41
fitting histogram for k = 51
fitting histogram for k = 61
fitting histogram for k = 71
table of predicted num. of genomic k-mers: histograms.dat
best k: 21

データ量が多いと計算時間もかなり大きくなるようで，出力にも少し書かれている通り，今回は1時間半くらいかかった．

さて，今回の場合はKmerGenieによる最適なk-merの値は21となった．出力結果で示されるk-mer頻度のプロットを並べて見てみると，kが大きい値の時には確かに正規分布として近似できないくらいなだらかな曲線になっているものの，K=21の最適値付近でもその様子はあまり変わっていないような感じがする．そもそも可視化の都合もあって横幅が揃えられていなかったりと色々と問題はあるが，とりあえずは今回のデータから推定された最適値は21だった．

KmerGenieで使用できるパラメータ

以上のように，KmerGenieは幾つかのk-merの値で実行して最適と思われる値を自動で出力するが，他にも手動でパラメータを指定することでk-merの範囲や刻みを変更することができる．

1
2
3
4
5
6
7
8
9
10
11
12
13
14

$ kmergenie
KmerGenie

Usage:
    kmergenie <read_file> [options]

Options:
    --diploid    use the diploid model
    --one-pass   skip the second pass
    -k <value>   largest k-mer size to consider (default: 121)
    -l <value>   smallest k-mer size to consider (default: 15)
    -s <value>   interval between consecurive kmer sizes (default: 10)
    -e <value>   k-mer sampling (default: auto-detected power of 10)"
    -o <prefix>  prefix of the output files (default: histograms)"

まとめ

このように，KmerGenieはシーケンスデータのみを用いてアセンブラに依存しない形で，アセンブルに指定する最適なk-merの値を求めることができる．推定された値が最良のアセンブルに繋がるかは実際に確認してみないと判断できないが，大まかな目安であったりk-mer決定プロセスの理由付けなどに力を発揮するだろう．今回は結果に対してそれほど詳細な評価を行わなかったが，arXiv.orgで公開されている論文の方ではGAGEのデータセットに対してアセンブル結果を含めた性能比較をしているので，そちらも参考にしていただきたい．

参考