トップページ > bash > Rで作成した行列データをいじりたおす
Rで作成した行列データをいじりたおす - bash
DNAマイクロアレイにしても次世代シークエンサーにしても、タブ区切りテキストをいじり倒すのがとても重要です
最近のデータはもはや、ECELで開けない量のデータが出てきています。そうでなくても、ECELでいちいちデータいじるのもしんどいですし。。
#Rの起動
$R
#各実験データ(wild type 1,2 vs mutant 1,2)のオブジェクトを生成する。
#オミックスデータの場合、要素数は万のオーダであることに注意。
wt1 <- rnorm(10,sd=5,mean=10)
wt2 <- rnorm(10,sd=5,mean=10)
mt1 <- rnorm(10,sd=5,mean=10)
mt2 <- rnorm(10,sd=5,mean=10)
#４つのオブジェクトを列方向で結合させて、ひとつのマトリックスオブジェクトを生成する。
data <- cbind(wt1,wt2,mt1,mt2)
#matrixからデータフレームに変換する。
data <- as.data.matrix(data)
#geneIDを格納したオブジェクトを生成する。
geneID <- c("gene1","gene2","gene3","gene4","gene5","gene6","gene7","gene8","gene9","gene10")
#geneIDとdataを結合させる。（geneIDすなわち行番号をcolnames(data)に組み入れず、意図的に列として扱うことで、シェルでも扱いやすくなる。)
data <- as.data.frame(geneID,data)
#表示してみる
data
geneID wt1 wt2 mt1 mt2
1 gene1 5.3540485 10.1280394 14.617105 6.858430
2 gene2 11.5216158 14.5030568 7.085197 13.102492
3 gene3 13.7710811 7.3251649 13.330784 10.927682
4 gene4 13.6855603 12.9628329 11.214203 8.208379
5 gene5 5.8908762 2.8704293 3.867760 13.401720
6 gene6 10.4802003 12.9800204 3.919811 13.068045
7 gene7 4.9378381 18.3811396 6.663597 8.449231
8 gene8 14.3298039 1.3903394 10.880709 10.668398
9 gene9 13.7565999 0.8564426 7.720071 7.457218
10 gene10 0.6470612 6.1940915 15.617938 18.897270
#タブ区切りテキストとして書き出す。row.names=Fであることに注意。
write.table(data, file="data.txt", append=F, quote=F, sep="\t",row.names=F, col.names=T)
#Rの終了
q()
#catコマンドでファイルの内容をざっと眺める。
$ cat data.txt
geneID wt1 wt2 mt1 mt2
gene1 5.35404848886385 10.1280394156775 14.6171048972252 6.85843001971764
gene2 11.5216157754471 14.5030567924109 7.08519743930237 13.1024919832259
gene3 13.7710811043585 7.32516494384991 13.330783789993 10.9276822183253
gene4 13.685560287049 12.9628329089167 11.2142032367318 8.20837940027492
gene5 5.8908761625754 2.87042929983391 3.86776035487622 13.4017204451707
gene6 10.4802002976654 12.9800203785410 3.91981136037804 13.0680452237444
gene7 4.93783813414918 18.3811395889717 6.66359654696842 8.44923140218957
gene8 14.3298038985442 1.39033944308932 10.8807085354603 10.6683978986145
gene9 13.7565999427345 0.856442626818601 7.72007073640271 7.4572175263714
gene10 0.647061248986092 6.1940915115274 15.6179381871434 18.8972700901370
#tailコマンドにより、ヘッダー（このファイルならば一行目）を飛ばして、２行目より最終行まで表示
$ tail -n +2 data.txt
gene1 5.35404848886385 10.1280394156775 14.6171048972252 6.85843001971764
gene2 11.5216157754471 14.5030567924109 7.08519743930237 13.1024919832259
gene3 13.7710811043585 7.32516494384991 13.330783789993 10.9276822183253
gene4 13.685560287049 12.9628329089167 11.2142032367318 8.20837940027492
gene5 5.8908761625754 2.87042929983391 3.86776035487622 13.4017204451707
gene6 10.4802002976654 12.9800203785410 3.91981136037804 13.0680452237444
gene7 4.93783813414918 18.3811395889717 6.66359654696842 8.44923140218957
gene8 14.3298038985442 1.39033944308932 10.8807085354603 10.6683978986145
gene9 13.7565999427345 0.856442626818601 7.72007073640271 7.4572175263714
gene10 0.647061248986092 6.1940915115274 15.6179381871434 18.8972700901370
#2列目(wt1)の合計値を計算してみる。
tail -n +2 data.txt | awk '{sum+=$2} END{print sum}'
94.3747
#2列目(wt1)の行数をカウントする。
$ tail -n +2 data.txt | awk '{ln++} END{print ln}'
10
#２列目(wt1)の平均値を求める。
$ tail -n +2 data.txt | awk '{sum+=$2 ; ln++} END{print sum/ln}'
9.43747
#各行の平均値を求める。
$ tail -n +2 data.txt | awk '{print ($2+$3+$4+$5)/4}'
9.23941
11.5531
11.3387
11.5177
6.5077
10.112
9.60795
9.31731
7.44758
10.3391
#各行においてFold Change(mt/wt)を求める(結果はリダイレクトによりファイルに書き出す。)
$ tail -n +2 data.txt | awk '{print ($4+$5)/($2+$3)}' > fc.txt
1.38712
0.775713
1.14989
0.728846
1.97111
0.724113
0.648091
1.3708
1.03861
5.04523
#Fold Changeの列用のヘッダーを記入したファイルを用意する。
echo 'Fold Change(mt/wt)' > header.txt
#結合(concatenate)させる。
$ cat header.txt fc.txt > fold_change.txt
$ less fold_change.txt
Fold Change(mt/wt)
1.38712
0.775713
1.14989
0.728846
1.97111
0.724113
0.648091
1.3708
1.03861
5.04523
#もとのデータファイル(data.txt)とFC計算値の格納したファイル(fold_change.txt)を横方向に結合(pasteコマンド）
$ paste data.txt fold_change.txt > output.txt
$ cat output.txt
geneID wt1 wt2 mt1 mt2 Fold Change(mt/wt)
gene1 5.35404848886385 10.1280394156775 14.6171048972252 6.85843001971764 1.38712
gene2 11.5216157754471 14.5030567924109 7.08519743930237 13.1024919832259 0.775713
gene3 13.7710811043585 7.32516494384991 13.330783789993 10.9276822183253 1.14989
gene4 13.685560287049 12.9628329089167 11.2142032367318 8.20837940027492 0.728846
gene5 5.8908761625754 2.87042929983391 3.86776035487622 13.4017204451707 1.97111
gene6 10.4802002976654 12.9800203785410 3.91981136037804 13.0680452237444 0.724113
gene7 4.93783813414918 18.3811395889717 6.66359654696842 8.44923140218957 0.648091
gene8 14.3298038985442 1.39033944308932 10.8807085354603 10.6683978986145 1.3708
gene9 13.7565999427345 0.856442626818601 7.72007073640271 7.4572175263714 1.03861
gene10 0.647061248986092 6.1940915115274 15.6179381871434 18.8972700901370 5.04523
#発現量の変化(Fod Change)の大きい遺伝子の順にソートをかける。
$ sort -k 6,6 -n -r output.txt
gene10 0.647061248986092 6.1940915115274 15.6179381871434 18.8972700901370 5.04523
gene5 5.8908761625754 2.87042929983391 3.86776035487622 13.4017204451707 1.97111
gene1 5.35404848886385 10.1280394156775 14.6171048972252 6.85843001971764 1.38712
gene8 14.3298038985442 1.39033944308932 10.8807085354603 10.6683978986145 1.3708
gene3 13.7710811043585 7.32516494384991 13.330783789993 10.9276822183253 1.14989
gene9 13.7565999427345 0.856442626818601 7.72007073640271 7.4572175263714 1.03861
gene2 11.5216157754471 14.5030567924109 7.08519743930237 13.1024919832259 0.775713
gene4 13.685560287049 12.9628329089167 11.2142032367318 8.20837940027492 0.728846
gene6 10.4802002976654 12.9800203785410 3.91981136037804 13.0680452237444 0.724113
gene7 4.93783813414918 18.3811395889717 6.66359654696842 8.44923140218957 0.648091
#発現量が低下(Fold Change<1)の遺伝子だけを抽出し表示する。
$ awk '$6<1 {print $0} output.txt
gene2 11.5216157754471 14.5030567924109 7.08519743930237 13.1024919832259 0.775713
gene4 13.685560287049 12.9628329089167 11.2142032367318 8.20837940027492 0.728846
gene6 10.4802002976654 12.9800203785410 3.91981136037804 13.0680452237444 0.724113
gene7 4.93783813414918 18.3811395889717 6.66359654696842 8.44923140218957 0.648091
#ちなみに、タブ区切りテキストであるこのファイルをRで読み込むには、read.delim()を用いる。
$ R
output <- read.delim("output.txt")
#outputオブジェクトの内容確認。
output
geneID wt1 wt2 mt1 mt2 Fold.Change.mt.wt.
1 gene1 5.3540485 10.1280394 14.617105 6.858430 1.387120
2 gene2 11.5216158 14.5030568 7.085197 13.102492 0.775713
3 gene3 13.7710811 7.3251649 13.330784 10.927682 1.149890
4 gene4 13.6855603 12.9628329 11.214203 8.208379 0.728846
5 gene5 5.8908762 2.8704293 3.867760 13.401720 1.971110
6 gene6 10.4802003 12.9800204 3.919811 13.068045 0.724113
7 gene7 4.9378381 18.3811396 6.663597 8.449231 0.648091
8 gene8 14.3298039 1.3903394 10.880709 10.668398 1.370800
9 gene9 13.7565999 0.8564426 7.720071 7.457218 1.038610
10 gene10 0.6470612 6.1940915 15.617938 18.897270 5.045230
