「ほっ」と。キャンペーン

<   2014年 01月 ( 2 )   > この月の画像一覧


2014年 01月 09日

makeblastdb

makeblastdb -in ref.fa -dbtype nucl -hash_index
blastn -db ref.fa -query que.fa -outfmt 6 > blastout.txt

makeblastdb -in ref.pep.fa -dbtype prot
blastp -db ref.pep.fa -query que.pep.fa -outfmt 6 > blastout.txt
[PR]

by koretoki | 2014-01-09 12:25
2014年 01月 07日

メタボデータをくっつけた

メタボローム解析の反復データをひとつにまとめたい。
n=3のものとn=1のものを比べたいのだけど、できればn=3全部使いたい。
(a1+a2+a3)/3 ってやっていいのかわからなかったから調べてみた。

#もともと持っているのはこんなデータ。いわゆるMartixとして得られたデータである。
name a1 a2 a3
CYT001 0.0007330321 0.0004886212 0.0008590029
CYT002 0.1763202000 0.1228157000 0.1730612000
CYT003 0.0137760700 0.0123575800 0.0164345100
CYT004 0.0121694300 0.0103101300 0.0095795930
CYT005 0.0008732506 0.0006902750 0.0008904565
CYT006 0.0007769990 0.0005277728 0.0005936377
CYT007 0.0010050660 0.0005964599 0.0010336850
CYT008 0.0000000000 0.0000000000 0.0000000000
CYT009 0.0000000000 0.0000000000 0.0000000000
CYT010 0.0005402510 0.0004378932 0.0003912269
CYT011 0.0000000000 0.0000000000 0.0000000000
CYT012 0.0000000000 0.0000000000 0.0000000000
#(これがずっと続いてる)

#R様に読んでいただく。
> data<-read.table("data.txt",header=T,row.names=1)
> attach(data)
> summary(data)
a1 a2 a3
Min. :0.0000000 Min. :0.0000000 Min. :0.0000000
1st Qu.:0.0003296 1st Qu.:0.0002691 1st Qu.:0.0000955
Median :0.0015751 Median :0.0015310 Median :0.0010649
Mean :0.0696770 Mean :0.0660012 Mean :0.0445601
3rd Qu.:0.0200691 3rd Qu.:0.0198329 3rd Qu.:0.0128484
Max. :2.4163290 Max. :2.2341410 Max. :1.5474730

#読んでいただけた。実際にはデータが全部表示されるが省略。
> data
a1 a2 a3
CYT001 0.0007330321 0.0004886212 0.0008590029
CYT002 0.1763202000 0.1228157000 0.1730612000
CYT003 0.0137760700 0.0123575800 0.0164345100
CYT004 0.0121694300 0.0103101300 0.0095795930
CYT005 0.0008732506 0.0006902750 0.0008904565
CYT006 0.0007769990 0.0005277728 0.0005936377
CYT007 0.0010050660 0.0005964599 0.0010336850

#総和がこんなにズレてる...
> sum(a1)
[1] 15.53798
> sum(a2)
[1] 14.71827
> sum(a3)
[1] 9.936896

#これをそのままプロットして1:1のところに線引いてみる。
> plot(a1,a3)
> abline(0,1,col="red")
e0160319_1135217.png

#ずれずれっぽいけどよく見えないからlog
> plot(log(a1),log(a3))
> abline(0,1,col="red")
e0160319_114586.png


#とりあえず総和で割ればいけそう
> plot(log(a1/sum(a1)),log(a3/sum(a3)))
> abline(0,1,col="red")
e0160319_1152369.png


(a1/sum(a1)+a2/sum(a2)+a3/sum(a3))/3
すればとりあえず物質間の比については使えるデータになりそう。
a1*x+a2*y+a3*zの形にしたい。
総和でも良さそうだけど怖いから回帰で比を出してもらうことにした。

#それぞれ回帰してもらう
#回帰モデルする
> fit_12<-lm(a1~a2)
> fit_23<-lm(a2~a3)
#中見る1vs2
> fit_12

Call:
lm(formula = a1 ~ a2)

Coefficients:
(Intercept) a2
-0.0005735 1.0643829
#中見る2vs3
> fit_23

Call:
lm(formula = a2 ~ a3)

Coefficients:
(Intercept) a3
0.006123 1.343758

#調節する時にかける値は
x:y=1:1.0643829
y:z=1:1.343758
x:y:z=1:1.0643829:1.430273=0.2861512:0.3045745:0.4092743

0.2861512*a1+0.3045745*a2+0.4092743*a3 で欲しい値が取れる。

#総和の比でやると
> sum(a1)
[1] 15.53798
> sum(a2)
[1] 14.71827
> sum(a3)
[1] 9.936896
> sum(a1)/sum(a2)
[1] 1.055693
> sum(a2)/sum(a3)
[1] 1.481174
> #x:y=1:1.055693
> #y:z=1:1.481174
> 1.055693*1.481174
[1] 1.563665
> #x:y:z=1:1.055693:1.563665
> 1*1/(1+1.055693+1.563665)
[1] 0.2762921
> 1*1.055693/(1+1.055693+1.563665)
[1] 0.2916796
> 1*1.563665/(1+1.055693+1.563665)
[1] 0.4320283
> #x:y:z=1:1.055693:1.563665=0.2762921:0.2916796:0.4320283

総和の方が強い?補正になってるのは、値の大きな子が引っ張ってるのかなー?

#回帰で求めた比がちゃんとハマるか心配
> plot(a1,a2)
> abline(0,1/1.0643829,col="red")
> plot(a2,a3)
> abline(0,1/1.343758,col="red")
> lm(a3~a1)

Call:
lm(formula = a3 ~ a1)

Coefficients:
(Intercept) a1
-0.00182 0.66564

> plot(a3,a1)
> abline(0,1/0.66564,col="red")

#logでもハマるか一応
> plot(log(a1),log(a2))
> e<-exp(1)
> e
[1] 2.718282
> abline(-log(1.0643829),1,col="red")
> plot(log(a2),log(a3))
> abline(-log(1.343758),1,col="red")
> abline(0,1,col="blue")
> plot(log(a3),log(a1))
> abline(0,1,col="blue")
> abline(-log(0.66564),1,col="red")

e0160319_116455.png


赤線が回帰の線
青線は1:1の線

人気ブログランキングに参加しています。
応援よろしくお願いします。
FC2 Blog Ranking

[PR]

by koretoki | 2014-01-07 01:18