2011-12-14 30 views
1

とクラスタデンドログラムの相関IはRで行われ、次のプロットを持っている: enter image description hereR:密度プロット

Iプロットを作るために、次のコードを使用する:

par(mfrow=c(1,2)) 

rmsd <- read.table(textConnection(" 
pdb rmsd 
1grl_edited.pdb 1.5118 
1oel_edited.pdb 1.1758 
1ss8_edited.pdb 0.8576 
1gr5_edited.pdb 1.8301 
1j4z_edited.pdb 0.7892 
1kp8.pdb 0.1808 
1kpo_edited.pdb 0.7879 
1mnf.pdb 1.2371 
1xck.pdb 1.6820 
2c7e_edited.pdb 5.4446 
2cgt_edited.pdb 9.9108 
2eu1.pdb 54.1764 
2nwc.pdb 1.6026 
2yey.pdb 61.4931 
"), header=TRUE) 

dat <- read.table(textConnection(" 
pdb  PA  EHSS 
1gr5_edited.pdb 21518.0 29320.0 
1grl_edited.pdb 21366.0 28778.0 
1j4z_edited.pdb 21713.0 29636.0 
1kp8.pdb 21598.0 29423.0 
1kpo_edited.pdb 21718.0 29643.0 
1mnf.pdb 21287.0 29035.0 
1oel_edited.pdb 21377.0 29054.0 
1ss8_edited.pdb 21543.0 29459.0 
1sx3.pdb 21651.0 29585.0 
1xck.pdb 21191.0 28857.0 
2c7e_edited.pdb 22930.0 31120.0 
2cgt_edited.pdb 22807.0 31058.0 
2eu1.pdb 22323.0 30569.0 
2nwc.pdb 21338.0 29326.0 
2yey.pdb 21032.0 28670.0 
"), header=TRUE, row.names=NULL) 

d <- dist(rmsd$rmsd, method = "euclidean") 
fit <- hclust(d, method="ward") 
plot(fit, labels=rmsd$pdb) 
groups <- cutree(fit, k=3) 

rect.hclust(fit, k=3, border="red") 

#for (i in dat[1]){for (z in i){ if (z=="1sx3.pdb"){print (z)}}} 

den.PA <- density(dat$PA) 
plot(den.PA) 
for (i in dat$PA){ 
    lineat = i 
    lineheight <- den.PA$y[which.min(abs(den.PA$x - lineat))] 
    lines(c(lineat, lineat), c(0, lineheight), col = "red") 
} 

左のプロットに示します右側のプロットは "PA"の密度プロットを示しています。参照がプロットに含まれていたので、密度プロットは、余分な値が含まれて明らかにそれはdatで参照ファイルは、クラスタのプロットがあり1sx3.pdb

で0の値を返しますので、参照はRMSDクラスタに含まれていませんでした3つの赤いボックスは、どうやって色を変えることができますか?左のボックスは赤、中央のボックスは緑、右のボックスは青です。私は密度プロットでミラーリングする必要があります。つまり、赤いボックスの中の値は密度プロット上に赤い線を持ち、緑のボックスの中の値は密度プロットなどの上に緑色の線を持っています。

参照構造をキャッチし、密度プロットで黒色にすることもできますか?

答えて

2

このコードは、あなたが望むことをします。あなたはほとんどそこにいました...ちょっとした並べ替えと索引付けが必要でした。

par(mfrow=c(1,2)) 

d <- dist(rmsd$rmsd, method = "euclidean") 
fit <- hclust(d, method="ward") 
plot(fit, labels=rmsd$pdb) 
groups <- cutree(fit, k=3) 

cols = c('red', 'green', 'blue') 

rect.hclust(fit, k=3, border=cols) 

#for (i in dat[1]){for (z in i){ if (z=="1sx3.pdb"){print (z)}}} 

cols = cols[sort(unique(groups[fit$order]), index=T)$ix] 

den.PA <- density(dat$PA) 
plot(den.PA) 
for (i in 1:length(dat$PA)){ 
    lineat = dat$PA[i] 
    lineheight <- den.PA$y[which.min(abs(den.PA$x - lineat))] 
    col = cols[groups[which(rmsd$pdb == as.character(dat[i, 'pdb']))]] 
    lines(c(lineat, lineat), c(0, lineheight), col = col) 
} 

enter image description here

+0

をシフトし、私は次のいずれかをトリミングしなければなりませんでしたデータセットが一致しなかったので、最初に 'dat'からの行。 –

+0

ちょっと、答えのためのthx。余分な行は、私の質問に記載されている参照構造です。 – Harpal

+0

作品は完璧に感謝:) – Harpal

0

あなたはこのようなように、国境に色のベクトルを渡すことができます。私はそれからの出力を保存し

t <- rect.hclust(fit, k=3, border=c("red",'green','blue')) 

注意、それは次のようになります。そして、

[[1]] 
[1] 12 14 

[[2]] 
[1] 1 2 3 4 5 6 7 8 9 13 

[[3]] 
[1] 10 11 

、あなたを変えますループを少しこれに

for (i in 1:length(dat$PA)){ 
    lineat = dat$PA[i] 
    lineheight <- den.PA$y[which.min(abs(den.PA$x - lineat))] 
    if(i %in% t[[1]]) lines(c(lineat, lineat), c(0, lineheight), col = "red") 
    if(i %in% t[[2]]) lines(c(lineat, lineat), c(0, lineheight), col = "green") 
    if(i %in% t[[3]]) lines(c(lineat, lineat), c(0, lineheight), col = "blue") 
} 

コードの最後のビットはあまりにもエレガントではありませんが、誰かがより良い解決策を考え出すことができると確信しています。

+0

プロットその後、上記の私の元のプロットはあなたが作る密度プロットをプロットした場合に密度のプロットは、一致していない、線の一部がちなみに – Harpal

関連する問題