case2 - nibb-unix/gitc202502-unix GitHub Wiki

実践演習２　RNA-Seq 解析結果の集計(2)

実践演習1 で作成したファイルは、各遺伝子について12 個のサンプルのリード数を記録したタブ区切りのテーブルであり、R で処理するのに適した形になっている。12 個のサンプルは、4つの条件について、それぞれ3 つずつのレプリケートをとったデータであった。

そこで、 R を用いて、各条件についてそれぞれ3 つのレプリケーションの平均値を計算してみよう。ただし、サンプルごとにリード数の総和が違っているため、それを補正するためにまず標準化（ノーマライズ）を行う必要がある。

ここでは、各サンプルのリード数の総和が100 万になるように標準化する CPM (Counts permillion reads) という値を計算しよう。それには、もとのリード数をサンプルごとに各サンプルのリード数の総和で割って100 万倍すればよい。なお、標準化については実践編で詳しく説明されるが、異なる遺伝子間で発現量を比較する場合は、リード数に加えて遺伝子の長さについての補正も必要で、そのために RPKM (Read per kilobase per million reads)という値がよく用いられる。これはもとのリード数を、遺伝子の長さ1,000 bp あたり、各サンプルのリード数の和100万あたりになるように標準化したものである。ここではより簡単なCPM を計算するが、同じ遺伝子の発現量をサンプル間で比較するだけの場合にはこれも有効な指標である。

以下でコマンド先頭の > はR のプロンプトを表すので、それ以降を改行なしで入力すること。

1.

まずR の作業ディレクトリを ~/デスクトップ/gitc/data/9_ex/rnaseq に移動しよう。メニュー（その他→作業ディレクトリの変更）から移動しても、コンソールから setwd("ディレクトリ名") を打ち込んでもよい。移動したら、getwd() で正しく移動できていることを確認しよう。

次に、ecoli.count_all ファイルからデータを読み込む。read.table 関数を使ってデータを読み込み、変数 eco_rna に代入しよう。このファイルは、セパレータはタブで、ヘッダはなしである。1 列目に遺伝子名が入っているので、これを各行の「名前」として指定しよう。これには、read.table のオプションとして row.names=1 を指定する。このとき、row.names で指定した列（1 列目）が行の名前として読み込まれる。

読み込んだら、head 関数で先頭数行を表示して、正しく読めているかを確認しよう。

最初の行で V1 がないことに注意。1 列目はデータではなく、行の名前として読み込まれている。行に名前をつけると eco_rna["b0020",] のように名前で行を抽出することが可能となり、興味のある遺伝子の情報を抽出したい場合などには便利である。

なお、1 行目はヘッダを表しており、最初のデータである b0001 は2 行目から始まる。「入力ファイルにヘッダあり」として読み込んでしまうと行がずれてしまうので注意しよう。また、今回はファイルにヘッダ行が含まれないので、各列にはデフォルトで R がつけた名前（V+元のファイルの列番号）がついている。列名を変更したい場合は colnames という関数で行えるが、ここでは省略する。

2.

CPM を計算するには、eco_rna の値をそれぞれのサンプルのリード数の和で割る必要がある。この和をまず計算しよう。eco_rna の各列に各サンプルのリード数が入っているので、それぞれの列を抽出したベクトルについて和をとるとよい。ベクトルの和をとるのは sum 関数で行えるので、たとえば eco_rna の1 列目の和は sum(eco_rna[,1]) で計算することができる（試して見よ）。この計算を各列について繰り返し行えばよい。

このように、行列やデータフレームから各行または各列をベクトルとして抽出して、それに特定の関数を適用する処理を繰り返すときは、apply 関数を使うのが便利である。apply は、apply(x, MARGIN, FUN) のように３つの引数をとる。x は入力データとなる行列またはデータフレーム、MARGIN は 1 または 2 をとり、1 は行（横）方向（図1A）、2 は列（縦）方向（図1B）のベクトルに分割した上で、行列全体にわたって処理を繰り返すことを指定する。FUN は適用する関数名である（help(apply) で確認せよ）。

この apply を使って、eco_rna の各列について列方向にベクトルを抽出し、sum 関数を適用して、結果を　eco_rna_readsum という変数に格納せよ。結果は要素数12 のベクトルになるはずである。その最初の要素が sum(eco_rna[,1]) と一致することを確認せよ。

3.

CPM を計算するには eco_rna の各行について、値を eco_rna_readsum の同じ位置の値で割って1,000,000 倍すればよい。今回は列方向ではなく、行方向にベクトルを抽出して計算することになるが、再び apply 関数を使うことができる。

まずここで行う計算を確認しよう。たとえば1 行目については以下のベクトルの割り算を行う。

> eco_rna[1,] / eco_rna_readsum

apply を使うには関数の形にする必要があるが、実は割り算の演算子 / は、クオートで囲むことによって2 つの引数を持つ関数名としても表すことができる。例えば、6 / 3 は、関数形で表すと、'/'(6,3) と書ける（試して見よ）。

従って，上記の計算を行うには apply に関数名として '/' を与えればよい。ただし、sum と違って '/' は引数を２つとる。１つめの引数は apply 内部で与えられるが、２番目の引数は外から与えなければならない。apply 関数では ４番目以降の引数に、適用する関数の２番目以降の引数を与えることができるので、除数 eco_rna_readsum を４番目の引数として指定する。こ れで行方向に関数を適用するように apply 関数を実行すればよい。apply を実行した後で全体を 1,000,000 倍する。この結果を eco_rna_cpm0 という変数に代入せよ。

これを head で確認すると、結果が大量に表示されて流れてしまうはずである。実は、ここでの結果は行と列が入れ替わってしまっており（これは apply 関数の仕様）、１行の長さが非常に長くなっている。これは dim 関数で確認できる。dim 関数は、指定した行列やデータフレームの行数と列数を表示する。dim(eco_rna) と dim(eco_rna_cpm0) を比較してみよう。

これでは扱いづらいので、eco_rna_cpm0 の行と列を入れ替えよう。これは転置行列をとる(transpose)という操作であり、R では t 関数で行う（help(t) で確認せよ）。この結果を eco_rna_cpm に代入しよう。head で結果を確認すると、以下のようになるはずである。

4.

最後に、標準化した値を使って、サンプルごとに３つのレプリケートの平均値を計算しよう。各行（遺伝子）について、各サンプルは、それぞれ1-3, 4-6, 7-9, 10-12 列目に対応しているので、それらの値の平均値を計算すればよい。

再び apply 関数を使おう。そのために、まず各行について、上記の４つの平均値を計算するための関数 calc_means を作る。関数の作成は、コマンドラインからそのまま打ち込んでもよいが、ちょっと複雑なのでエディタを使って行うとよいだろう。

calc_means <- function(v) {
  c( mean(v[1:3]), mean(v[4:6]), mean(v[7:9]), mean(v[10:12]) )
}

この関数は、引数に v として各行における要素１２のベクトル値をとり、それを３要素づつの４つのベクトルに分けて、それぞれの平均値を計算し、c() 関数によって値をまとめたベクトルとして返している。calc_means(1:12) を実行して動作を確認しよう。

前問と同様にして、apply 関数によってこの関数を eco_rna_cpm に適用し、サンプルごとの平均値を計算しよう。その際、行と列の入れ換えにも注意すること。正しく作成できれば、head で確認すると、以下のような結果が得られるはずである。