BioConductorのGenomicFeaturesパッケージにあるTranscriotDbオブジェクトを使う際に関連したクラスとの相互関係がよく分からなくなってきたので，まとめました．

> sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
locale:
[1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] BSgenome.Hsapiens.UCSC.hg19_1.3.19 BSgenome_1.26.1
[3] Biostrings_2.26.3 TxDb.Hsapiens.UCSC.hg19.knownGene_2.8.0
[5] GenomicFeatures_1.10.2 AnnotationDbi_1.20.7
[7] Biobase_2.18.0 GenomicRanges_1.10.7
[9] IRanges_1.16.6 BiocGenerics_0.4.0
loaded via a namespace (and not attached):
[1] biomaRt_2.14.0 bitops_1.0-4.2 DBI_0.2-5 parallel_2.15.1
[5] RCurl_1.95-4.1 Rsamtools_1.10.2 RSQLite_0.11.2 rtracklayer_1.18.2
[9] stats4_2.15.1 tools_2.15.1 XML_3.95-0.2 zlibbioc_1.4.0

2012-12-14

Perl でawkライクなワンライナー（GFFからの行抽出など）

Perl

メモ．awkでいいという噂もありますが…

参考

2012-12-12

RStudioでknit HTMLで図が出力されなかったのを解決した

いつの間にか，RStudioの knit HTMLボタンで作ったHTMLに図が出力されなくなったので（問題の分析が出来てないまま@wakutekaに相談したり…すみません…）色々試したり調べたりしていました．

こちらの掲示板で，knitrはevaluateに依存しているけど，そのevaluateだけが新しいバージョンになっていると，図が出力されなくなる（？）からknitrもアップデートしなよ，と書いてありました．
knitr - Knit does not produce figures - Stack Overflow

sessionInfo() で確かめると以下の組み合わせになっていました．

knitr_0.8
evaluate_0.4.3

install.packages("knitr") でインストールし直すと，以下のようになり，図の出力もされるようになりました．

knitr_0.9
evaluate_0.4.3

2012-12-10

Bowtie2でmultiple hitの扱いがおかしい件（Tophat，Cufflinksとの関連で）

tool alignment ngs RNA-seq

gaou_akさんのツイートで知りました．

ATTENTION: bowtie2 and multiple hits - SEQanswers

「Bowtie2のマニュアルでは，best hitが複数あるとランダムに選ぶ」と書いてあるにも関わらず，デフォルトのパラメタでは全てのリードが同じ位置にマップされるそうです（これはリアルに困る…）．
対策としては，すべてのmultiple hitを出力させるオプション（-k か -a）を使うということのようです．

multi-hit readsを除くべきか否か

上の掲示板にあったコメントで，「全て除く」「全て残す」「残して重み付け」の3つの方法が紹介されていました，

In my opinion, you can do 3 different things for downstream analysis:
1. Discard all reads with multiple hits
For a DE analysis it should be ok, since you delete them for all runs. But you will underestimate, or even completely miss the expression of all genes with multiple copies.
2. Take all hits
For a DE analysis it should be ok, since you do the same for all runs and you won't miss genes with multiple copies. But you will overestimate the expression of all genes with multiple copies.
3.Take all hits, but divide the "expression impact" of the read by the number of its mapping loci
I think this is the prefered way, since you do not know where the read actually comes from. You have all mappings, but by dividing the impact, you will "punish" the hit and equally distribute its impact in expression (it's somehow the same result as the random generator, but in my opinion much cleaner). That is actually, how we handle this issue by default. This way, you minimize the error and don't loose any region, where the read might come from. Another benefit of this approach is the possibility that you can always access the mapping loci at a later time and filter out multiple hits, or change the way of how you want to normalize them. Of course, most of the downstream tools are not able to handle this normalization, but we use our own tools for that.
ATTENTION: bowtie2 and multiple hits - SEQanswers

Cufflinksのマニュアルによると，Cufflinksは上記の3. の方針のようです．

デフォルトだと，アラインメントの数で割った寄与に重み付けされる（上記の3.）
さらに -u/--multi-read-correct を使うと，発現量推定の結果に合わせて重み付けが更新される

Tophat2のデフォルトではどうなっているか

Tophat2 2.0.6 でのログ（logs/run.log）を確認した所，デフォルトで動かした場合，bowtie2を「-k 」オプションで動かすようです．
また，tophat2のマニュアルによると，デフォルトでは「-g/--max-multihits」が20で，「--report-secondary-alignments」がOFFになっているため，20箇所まではbest alignment scoreをもったアラインメントが出力されます．
SAMのFLAGでは，best alignment scoreを持つアラインメントのうち，１つ以外は 0x100（secondary alignment）が立つそうですが，alignment scoreとしてはprimaryなものと同一だと考えられます．

2012-12-07

確率変数としてのp-valueの感覚をRで身につける

R statistics

R Advent Calendar 2012の7日目です．もう一週間経ったのですね．

The American Statistician, August 2008, Vol. 62, No. 3 P-Values are Random Variables (PDF) で「学生がp-valueについて理解するためにはp-valueが確率変数であることを強調した方がいいんでない」という提案がありました．
そこで，何の脈絡もありませんが，p-valueは確率変数であり，帰無分布からのp-valueは一様分布に従うことをRで簡単に確かめたいと思います．

帰無分布からのp-valueは一様分布に従う

まず，ある統計検定量Tが帰無仮説H0の下で確率密度関数 f(T) に従うとします．
p-valueの定義は，帰無仮説H0の下である検定統計量Tがt以上の値を取る確率です．
$p=P(T \gt = t|H_{0})$ (1)
p-valueは以下のようにも書けます．
$p=1-P(T \lt t|H_{0})$ (2)
ここで，累積分布関数 F(T) を考えると，次のように書けます．
$p=1-F(t)$ (3)
F(T)は単調増加なので，
$P(T \lt t|H_{0})=P(F(T) \lt F(t))$ (4)
(2), (3), (4) より
$P(F(T) \lt F(t)) = F(t)$ (5)
(5) より，F(t)は一様分布となり，1-F(t)も一様分布となるので，p-valueも一様分布となります．

Rで実際に確かめる

（確率に関連する関数について詳しくはこちらをご覧ください scratch-R: probability plots）

まず，ある統計量が平均1，分散1の正規分布に従うとします（帰無仮説）．そこから10000個のデータをサンプリングします．これはある統計検定（棄却されなそうな）を10000回やっているものと考えて下さい．

m <- 1  #平均
v <- 1  #分散
n <- 10000  #データ数
data <- rnorm(n, m, v)  #乱数生成

次に，この統計量が平均1，分散1の正規分布に従うときのp-valueの分布を描くと，一様分布になります．

pvalues <- 1 - pnorm(data, m, v)
hist(pvalues)

では，「平均0.5，分散1の正規分布」もしくは「平均2，分散1の正規分布」に従うときのp-valueの分布はどうなるでしょうか．

m.fake1 <- 0.5
pvalues.fake1 <- 1 - pnorm(data, m.fake1, v)
hist(pvalues.fake1)

m.fake2 <- 2
pvalues.fake2 <- 1 - pnorm(data, m.fake2, v)
hist(pvalues.fake2)

偏った分布になります．これは，さまざまな対立仮説ではp-valueの値のとりやすさが変わるということを直感的に表しています．

参考

2012-12-03

BAMヘッダーのサイズが大きいとCufflinksでエラーが出る

RNA-seq tool ngs

このリンク先のSEQanswersの掲示板は2011年下旬のものなのですが，Cufflinks v2.0.1, v2.0.2 でも同様の現象を確認したので，ここにまとめます．

エラー

Cufflinksを実行したら，以下のようなエラーが出ました．

You are using Cufflinks v2.0.2, which is the most recent release.
Warning: BAM header too large
File ./accepted_hits.bam doesn't appear to be a valid BAM file, trying SAM...
[16:47:02] Inspecting reads and determining fragment length distribution.
SAM error on line 175: CIGAR op has zero length
SAM error on line 177: CIGAR op has zero length
SAM error on line 196: CIGAR op has zero length
SAM error on line 218: CIGAR op has zero length
SAM error on line 220: CIGAR op has zero length
SAM error on line 267: CIGAR op has zero length
...

原因

Cufflinks のhits.cppにある以下の行によって，BAMのヘッダーが4MB以上の場合，BAMのパースができなくなるのだそうです．

static const unsigned MAX_HEADER_LEN = 4 * 1024 * 1024;
Cufflinks, BAM header problem solved... for the moment - SEQanswers

BAMのパースに失敗するとCufflinksは入力ファイルをSAMと判断して読もうとしますが，実際にはBAMなので，当然エラーが出て実行が停止します．

実際に自分のBAMのヘッダーのサイズを確認した所，11MBもありました．

samtools view -H hoge.bam > header.txt
ls -lsh header.txt

対策

hits.cppの数字を書き換えて再コンパイル
BAMをSAMに変換してCufflinksの入力ファイルとする

今回は面倒だったので，後者の方法をとることで解決しました．

samtools view -h hoge.bam > hoge.sam

教訓

ヒトやマウスを扱うことが多かったので，まさかBAMのヘッダーが4MBを越えているとは思いませんでした．
しかし，非モデル生物のドラフトゲノムではscafoldが1万から100万にもなるケースもあるので，今後このエラーにぶつかる方も増えると予想されます．

2011-12-08

ゲノム座標の記法（1-based vs. 0-based）

format

先日、後輩から「UCSCからダウンロードした遺伝子の座標から塩基配列を取得するときに、ファイルに書かれた座標とゲノムブラウザでみた配列が1つずれてる」と言われました。
自分もかつて同じ問題にぶつかったのですが、意外と知らない人がいると思い、まとめました。

まず、座標の表現方法には 0-based と1-based があります。0-basedは配列の一番最初の座標を0として扱います。1-basedは配列の一番最初の座標を1とします。一般にバイオ系の人は1-based、インフォ系の人は0-basedをよく使う気がします。

UCSCゲノムブラウザでよくみる「chr1:111111-222222」などという表示は1-basedです。一方、プログラムで座標を扱う場合には、0-basedが便利だったりします。

実はUCSCのFAQにはこんなことが書いてあります（意訳）。

ウチのデータベースの内部表現は、startは0-based、endは1-basedなんだよ。ゲノムブラウザで表示するときは、1足しているんだよ。
なんで、そんなめんどくさいことするかって？そりゃあ、単純な計算で長さとかが計算できた方が便利だからだよ。
http://genome.ucsc.edu/FAQ/FAQtracks#tracks1

例えば、「chr1:111111-222222」という領域を表したい場合、UCSCのデータベースからダウンロードしたファイルではこのようになります。

chr1 111110 222222

なので、もしこの領域の塩基配列を取得したい時には、startに1を足して扱わなければなりません。

これは、遺伝子のstrandが(-)の場合でも変わりません。つまり、Watson strandからみたときのstartとendとして表示されているので、startが0-basedとなります。

chr1 111110 222222 +
chr1 111110 222222 -

これらは同じ領域（chr1:111111-222222）のWatson strand側とCrick strand側を表します。

BEDでも、0-based start & 1-based endになっています。GTFやGFFは、1-based start & 1-based endです。
種々のフォーマットにおける座標の扱いについてもUCSCのFAQに載ってます( http://genome.ucsc.edu/FAQ/FAQformat.html )。

micro*macro

TranscriotDbオブジェクトと周辺のクラスとの関係