【R】遺伝子発現解析を学んでみたい②

・お題：前回の続き。Differential expression analysis（DEG）という解析を少し齧ってみたい。私はWet専門で学んできた人間だから、この辺の解析方法にはとんと疎い。以下のHPを順になぞって学んでいきたいと思う。

combine-australia.github.io

・今回は前回の続きであり、データセットyができている状態とする。参考サイトでいうと、Filtering lowly expressed genesの直前まで終わっている感じ。

> y

・次に、いずれのサンプルでも発現量が低すぎる遺伝子を除く。発現量が低すぎる遺伝子の除き方はいくつかあるらしいけれど、今回はcpmを計算して0.5未満のものを除外しているらしい。では、cpmって何か？cpmはedgeRのcpm関数で算出されるので、公式サイトを除いてみる。

rdrr.io

・Compute counts per million (CPM)とのこと。もっと言うと、サンプルごとに補正した各遺伝子のカウント数、という感じ。

・例えば、MCL1.DGを例にとると、lib.sizeが23227641なので、読めた配列が23227641個だったことになる（重複あり。固有の配列という意味ではない）。サンプルごとに状態が異なることなどから、読めるリード数も異なる値になる。単純に各遺伝子に関して読まれた数だけで比較すると、サンプルの状態のバイアスがかかってしまうので、読めた数で割って補正をかけてやるのが一般的らしい。

・読めた総数は数百万～数千万に及ぶことが一般的なのに対して、各遺伝子ごとのカウントは非常に小さいことになるので、million（百万）倍することで、見やすくするっぽい。

・ということで、とりあえず、やってみる。

> myCPM <- cpm(countdata)

> View(myCPM)

・cpmの閾値を0.5に設定し、各サンプルごとに何個の遺伝子が閾値を超えたか見てみる。

> thresh <- myCPM > 0.5

> table(rowSums(thresh))

・これで、TRUEが0個が10857遺伝子、1個が518遺伝子、2個が544遺伝子、、、12サンプルすべてでTRUEになっているのが11433遺伝子あったことになる。

・ここで、cpmの閾値を0.5に設定した理由が気になる。これは、カウントが10～15を超えた遺伝子、という意味みたい。ライブラリが大きくなると、閾値を小さくするらしい。

・少なくとも2つのサンプルで0.5以上だった遺伝子をキープしたい。

> keep <- rowSums(thresh) >= 2

> head(keep)#keepの最初の方を覗いてみる。

> summary(keep)#条件を満たす遺伝子がいくつあったか確認する。先の544から11433を全部足したのと同じTRUEの数になる。
Mode FALSE TRUE
logical 11375 15804

・これで解析対象にする遺伝子群のフィルターができたことになる。yにフィルターを適用し、解析対象のデータを抜き出す。

> y <- y[keep, keep.lib.sizes=FALSE]

> y#yを覗いてみる。

・確かに、yを作成した際にあったいずれのサンプルでも発現量が少ない遺伝子がなくなり、どれかのサンプルでそれなりに発現のある遺伝子ばかりになっている。

・次に、このデータセットのクオリティを見てみる。

・まずは、Libararyサイズを見てみる。それぞれのサンプルでどれだけの数の遺伝子（重複アリ）を読み取れたかを確認する。リード数の合計を見ればよく、yのsamplesのlib.sizeに格納されている。サンプルごとのリード数を棒グラフで見てみる。せっかくなので、ggplot2の練習も。

> library(tidyverse)

> ggplot(data = y$samples,
+ aes(x=row.names(y$samples),
+ y=y$samples$lib.size/1e06,
+ fill = row.names(y$samples)))+
+ geom_bar(stat = "identity")+
+ labs(title = "Barplot of library sizes",
+ x = "Samples",
+ y = "Library size (millions)")

・ぱっと見それほど顕著な差はない。なんとなく塗ってみたが、却ってわかりにくい…。。次に、各サンプルの遺伝子の発現量の分布をみてみる。遺伝子の発現量は、先にcpmとしてサンプルごとに正規化されたものがベースになるが、分布が正規分布ではなく、log2して対数変換してから図示する。

> logcounts <- cpm(y,log=TRUE)

> ggplot(data = logcounts %>% as.data.frame() %>% stack(),
+ aes(x = ind,y = values))+
+ geom_boxplot()+
+ geom_hline(aes(yintercept = median(logcounts)), color = "red")+
+ labs(title = "Boxplots of logCPMs (unnormalised)",
+ x = "Samples",
+ y = "Log2 counts per million")