Rのバッチ処理で三角形の面積を計算させてみた。

$ vim triangle.R

takasa <- 30
teihen <- 4
menseki <- takasa * teihen /2
menseki


$ R CMD BATCH triangle.R triangle_R.txt



$ cat triangle_R.txt

R version 2.10.1 (2009-12-14)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

Rは、自由なソフトウェアであり、「完全に無保証」です。
一定の条件に従えば、自由にこれを再配布することができます。
配布条件の詳細に関しては、'license()'あるいは'licence()'と入力してください。

Rは多くの貢献者による共同プロジェクトです。
詳しくは'contributors()'と入力してください。
また、RやRのパッケージを出版物で引用する際の形式については
'citation()'と入力してください。

'demo()'と入力すればデモをみることができます。
'help()'とすればオンラインヘルプが出ます。
'help.start()'でHTMLブラウザによるヘルプがみられます。
'q()'と入力すればRを終了します。

> takasa <- 30
> teihen <- 4
> menseki <- takasa * teihen /2
> menseki
[1] 60
>
>
> proc.time()
   ユーザ   システム       経過 
     0.460      0.052      1.309

Cで変数の格納

$ vim practice2.c

a


#include <stdio.h>

int main(void)
{
int suti1, suti2, suti3;
suti1 = 10;
suti2 = 20;
suti3 = 30;
printf("suti1の値は%dです。suti2値は%dです。suti3の値は%dです\n", suti1, suti2, suti3);

return 0;
}

Esc

:wq

$ gcc -o practice2 practice2.c
$ ./practice2
suti1の値は10です。suti2値は20です。suti3の値は30です


int型変数suti1,suti2,suti3を宣言してから、各々に整数を格納し、print関数によって表示させた。

私の専攻はbioinformaticsです。(C language)

$ vim practice1.c



a

#include<stdio.h>

int main(void)
{
printf("私の専攻は\n");
printf("Bioinformaticsです\n");
return 0;
}


Esc

:wq


$ gcc -o practice1 practice1.c
$ ./practice1
私の専攻は
Bioinformaticsです

正規分布に従う乱数の作成及びプロット。シミュレーション。

$ R

install.packages("scatterplot3d")

library(scatterplot3d)

a <- rnorm(10) # 正規分布に従う乱数を10個生成
b <- rnorm(10) # 正規分布に従う乱数を10個生成
c <- rnorm(10) # 正規分布に従う乱数を10個生成

d <- rnorm(100) # 正規分布に従う乱数を100個生成
e <- rnorm(100) # 正規分布に従う乱数を100個生成
f <- rnorm(100) # 正規分布に従う乱数を100個生成

g <- rnorm(1000) # 正規分布に従う乱数を1000個生成
h <- rnorm(1000) # 正規分布に従う乱数を1000個生成
i <- rnorm(1000) # 正規分布に従う乱数を1000個生成

x <- rnorm(10000) # 正規分布に従う乱数を10000個生成
y <- rnorm(10000) # 正規分布に従う乱数を10000個生成
z <- rnorm(10000) # 正規分布に従う乱数を10000個生成

png("110217.rnorm.png")
par(mfrow = c(2,2))
scatterplot3d(a,b,c)
scatterplot3d(d,e,f)
scatterplot3d(g,h,i)
scatterplot3d(x,y,z)
dev.off()


う~ん。胸が熱くなる。

scatterplod3dで3次元散布図を作成した★

$ R

install.packages("scatterplot3d")

library(scatterplot3d)


x <- c(1,2,3,4,5,6,7,8,9)
y <- c(9,8,7,6,5,4,3,2,1)
z <- c(21,22,23,24,25,26,27,28,29)

png("110217.scatterplot3d.png")
scatterplot3d(x,y,z)
dev.off()


perlで塩基配列を扱ってみたよん(^^)

perlを用いて、DNAの塩基配列を出力したり、二つのDNAフラグメントを結合したりしてみました(^_^)

$ vim 110215.1.pl

a

use strict;
use warnings;
#Storing DNA in a variable, and printing it

#First we store the DNA in a variable called $DNA
my $DNA = 'AGCGGAGGACGGGAAAATTACTACGGCATTAGC';

# Next, we print the DNA onto the screen
print "$DNA\n";



Esc

:wq



$ perl 110215.1.pl
AGCGGAGGACGGGAAAATTACTACGGCATTAGC



$ vim 110215.2.pl

a

  1 use strict;
  2
  3 #Store two DNA fragments into two vairiables called $DAN1 and $DNA2
  4 my $DNA1 = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC';
  5 my $DNA2 = 'ATAGTGCCGTGAGAGTGATGTAGTA';
  6
  7 # Print the DNA onto the screen
  8 print "Here are the original two DNA fragments:\n\n";
  9 print "DNA1 is $DNA1\n";
 10 print "DNA2 is $DNA2\n\n";
 11
 12 #Concatenate the DNA fragments into a third variable and print them
 13 #Using #string interpolation"
 14 my $DNA3 = "$DNA1$DNA2";
 15
 16 print "Here is the concatenation of the first two fragments(version1):\n\n";
 17 print "$DNA3\n\n";
 18
 19 #An alternative way using the "dot operator";
 20 #Concatenate the DNA fragments into a third variable and print them
 21 my $DNA3 = $DNA1 . $DNA2;
 22
 23 print "Here is the concatenation of the first two fragments(version2):\n";
 24 print "$DNA3\n\n";
 25
 26 #print the same thing without using the variable $DNA3
 27 print "Here is the concatenation of the first two fragments(version3):\n";
 28 print $DNA1, $DNA2, "\n";
 29
 30


Esc

$ perl 110215.2.pl
Here are the original two DNA fragments:

DNA1 is ACGGGAGGACGGGAAAATTACTACGGCATTAGC
DNA2 is ATAGTGCCGTGAGAGTGATGTAGTA

Here is the concatenation of the first two fragments(version1):

ACGGGAGGACGGGAAAATTACTACGGCATTAGCATAGTGCCGTGAGAGTGATGTAGTA

Here is the concatenation of the first two fragments(version2):
ACGGGAGGACGGGAAAATTACTACGGCATTAGCATAGTGCCGTGAGAGTGATGTAGTA

Here is the concatenation of the first two fragments(version3):
ACGGGAGGACGGGAAAATTACTACGGCATTAGCATAGTGCCGTGAGAGTGATGTAGTA

BMI, MCV, MCH, MCHCを計算&正常判定をRでやってみた。

患者さんの検査値が以下のようだとします。

Wight 160cm
Height 47kg

WBC 6800/ul
RBC 412*10^4/ul
Hct 38.8%
Hb  11.8g/dl


$ R

Height <- 160*10^-2
Weight <- 47

BMI <- Weight / Height^2      #Body Mass index

BMI
[1] 18.35937   #22ぐらいが丁度いいと言われている。




WBC <- 6800     #white blood cell.白血球。
RBC <- 4.12 * 10^6  # Red blood cell。赤血球。
Hct <- 38.8  #Hct:ヘマトクリット。全血二体する赤血球容積。
Hb  <- 11.8 #ヘモグロビンの濃度。

(4300<WBC) && (WBC<6900)
[1] TRUE

(446*10^4<RBC) && (RBC<515*10^4)
[1] FALSE


(43.2<Hct) && (Hct<48.6)
[1] FALSE


MCV <- Hct*10*10^6/RBC    #平均赤血球容積

 MCV
[1] 94.17476


(91.3<MCV) && (MCV<99.3)
[1] TRUE


MCH <- Hb*10*10^6/RBC #平均赤血球ヘモグロビン量

MCH
[1] 28.64078


(30.0<MCH) && (MCH<32.7)
[1] FALSE


MCHC <- Hb*100/Hct   #平均赤血球ヘモグロビン濃度

MCHC
[1] 30.41237

(32.3<MCHC) && (MCHC<33.7)
[1] FALSE

Windows XPで未割り当て領域にCドライブを拡張してみた!!

” EASEUS Partition Manager Home Edit  ”

は結構すごいですよ!!

なんといっても、マウスを使ってグニャッと動かすだけでパーティションの両々を変更できちゃう!!

http://www.partition-tool.com/personal.htm

にいって









解凍は、Lhasaみたいな解凍ソフトで一瞬で解凍できます。

中にインストーラー入ってるんで、さくさくインストールします。
インストールしたら、実際に


 EASEUS Partition Manager Home Editを立ち上げてみます。



赤い波線で囲ったパーティションの間のエッジをグイグイ動かせばパーティションの容量を自由自在に変更することができます♪






最後は、applyをクリックすれば変更が反映されました★☆

Windows環境のRで環境変数をいじってRgraphvizのインストールに成功したよ!!♪

http://www.graphviz.org/pub/graphviz/stable/windows/






ここで

graphviz-2.20.3a.msi

をダウンロードしてインストール!!!


次に、
コントロールパネル -> システム -> 詳細設定タブ -> 環境変数(N)

すると以下のウィンドウが開く。







ここで、下段のシステム環境変数(S):赤線でかこってあるやつ:をクリック。


変数名  Path 
変数値   C:\Program Files\Graphviz2.20\bin\

を入力。

あとは、2回分OKを押す♪



一度、パソコンを再起動かける♪再起動しないとうまくいかなかったわ(涙)



これで、Windows XPに対する設定は終わり!

最後にRを開いて、パッケージ Rgraphvizをインストールしてみる☆


source("http://bioconductor.org/biocLite.R")
biocLite("Rgraphviz") 



library("Rgraphviz")


 要求されたパッケージ Biobase をロード中です

Welcome to Bioconductor

  Vignettes contain introductory material. To view, type
  'openVignette()'. To cite Bioconductor, see
  'citation("Biobase")' and for packages 'citation(pkgname)'.

 要求されたパッケージ Category をロード中です
 要求されたパッケージ AnnotationDbi をロード中です
 要求されたパッケージ DBI をロード中です







# キタ━━━(゚(゚∀(゚∀゚(☆∀☆)゚∀゚)∀゚)゚)━━━!!!


#動作テストしてみます。

set.seed(123)
V  <-  letters[1:10]
M  <-  1:4
graph  <-  randomGraph(V,  M,  0.2)
graph  <-  layoutGraph(graph)

png("110214.graph.png")
renderGraph(graph)
dev.off()




Gene ontologyのデータをGO.dbで取得したのちネットワークを描かせた♪♪

GO, itself, is a structured terminology. The ontology describes genes and gene products and is divided into three separate ontologies. One for cellular component (CC),one for molecular function (MF) and one for biological process (BP).

つまり、

BPは、生物学的プロセス
MFは、分子機能
CCは、細胞の構成単位

っていうことみたいね。

これらのGOのデータベースを覗いて見ようと思う。

install.packages("igraph")
source("http://bioconductor.org/biocLite.R")
biocLite("GO.db")


library(GO.db)
library(igraph)


BP <- toTable(GOBPPARENTS)
CC <- toTable(GOCCPARENTS)
MF <- toTable(GOMFPARENTS)

#これで、BP(biological process)と、CC(Celular Component)と、MF(Molecular Function)をデータテーブルにできた!!!

#そしたら、これらを覗いてみよう♪

BP[1,]
       go_id      go_id RelationshipType
1 GO:0000001 GO:0048308              isa
CC[1,]
       go_id      go_id RelationshipType
1 GO:0000015 GO:0043234              isa
MF[1,]
       go_id      go_id RelationshipType
1 GO:0051082 GO:0005515              isa

#たとえば、BPの1から1000行までの一列目と二列目を抜き出してオブジェクトBP2を生成すると

BP2 <- BP[1:1000, 1:2]
BP2[1,]
       go_id    go_id.1
1 GO:0000001 GO:0048308

class(BP2)
[1] "data.frame"

BP3 <- as.matrix(BP2)


#ここで、BP3オブジェクトは、辺リストとみなすことができるので、graph.edgelist(辺リストの行列)により、igraphで利用可能なgraphオブジェクトを生成することができる♪


g.BP3 <- graph.edgelist(BP3)

plot(g.BP3, vertex.size = 1, vertex.color = 5, edge.color = 7, edge.arrow.size = 0.1 vertex.label.cex = 0.001, )


#vertex.size :頂点の大きさ。デフォルトは15。vertex.color:頂点の色。vertex.label.cex:頂点のラベルの大きさ.edge.arrow.size:矢印の大きさ。デフォルトは1。

png("test.png")
plot(g.BP3, vertex.size = 3, vertex.color = "yellow", edge.color = "green", edge.arrow.size = 0.01, vertex.label.cex = 0.1)
dev.off()




igraphでランダムグラフの作成♪♪

ランダムグラフとは、グラフに含まれる2つの頂点の間に一定の確率で辺を張ることによって作られるグラフである。ランダムグラフについて最初に体系的な研究をまとめた2人の数学者の名前から、Erdos-Renyi graphとも呼ばれる。

random.graph.game()の引数は

頂点数
2頂点間で辺を張る確率p
無向グラフ(デフォルト)か有向グラフかを決める。


例えば、

$ R

library(igraph)

g1 <- random.graph.game(50, p=0.01)
g2 <- random.graph.game(50, p=0.01, directed = TRUE)
g3 <- random.graph.game(50, p=0.1)
g4 <- random.graph.game(50, p=0.1, directed = TRUE)


png("110214.random.graph.game.png")
par(mfrow = c(2,2))
plot(g1)
plot(g2)
plot(g3)
plot(g4)
dev.off()



igraphで無向グラフを描写する。

#無向グラフを作成する

e.matrix <- matrix(c(1,2,1,3,1,4,2,3),ncol =2 , byrow = TRUE)


e.matrix
     [,1] [,2]
[1,]    1    2
[2,]    1    3
[3,]    1    4
[4,]    2    3

g4 <- graph.edgelist(e.matrix-1, directed = FALSE)

g4
Vertices: 4
Edges: 4
Directed: FALSE
Edges:
         
[0] 0 -- 1
[1] 0 -- 2
[2] 0 -- 3
[3] 1 -- 2

png("110214.g4.png")
plot(g4)
dev.off()



辺リストを与えて有向グラフを描く♪

#ベクトルとして辺リストを作成する。
e.list <- c(1,2,1,3,2,1,2,3,4,1)




#行列として辺リストを作成する。
e.matrix <- matrix(e.list, ncol = 2, byrow = TRUE)

#このbyrow=TRUEのテクニックは必見・・・。普通かな笑

e.list
 [1] 1 2 1 3 2 1 2 3 4 1

e.matrix
     [,1] [,2]
[1,]    1    2
[2,]    1    3
[3,]    2    1
[4,]    2    3
[5,]    4    1

e.list-1
 [1] 0 1 0 2 1 0 1 2 3 0


g2 <- graph(e.list-1, n=4)


e.matrix-1
     [,1] [,2]
[1,]    0    1
[2,]    0    2
[3,]    1    0
[4,]    1    2
[5,]    3    0



g3 <- graph.edgelist(e.matrix - 1)


png("110214.g23.png")
par(mfrow = c(1,2))
plot(g2)
plot(g3)
dev.off()


igraphで隣接行列を与えて有向グラフを作成してみた!!

$ R


library(igraph)



dg1 <- matrix(c(0,1,1,0,1,0,1,0,0,0,0,0,1,0,0,0), nrow=4, ncol = 4, byrow = TRUE)

dg1
     [,1] [,2] [,3] [,4]
[1,]    0    1    1    0
[2,]    1    0    1    0
[3,]    0    0    0    0
[4,]    1    0    0    0

class(dg1)
[1] "matrix"

colnames(dg1) <- 1:4
rownames(dg1) <- 1:4

dg1
  1 2 3 4
1 0 1 1 0
2 1 0 1 0
3 0 0 0 0
4 1 0 0 0


#igraphのパッケージでは、隣接行列をそのままネットワークデーとして用いることができないので、隣接行列からigraphで利用可能なigraphオブジェクトを生成する必要がある。


g1 <- graph.adjacency(dg1)

#adjacencyとは”隣接”という意味!


g1
Vertices: 4
Edges: 5
Directed: TRUE
Edges:
             
[0] '1' -> '2'
[1] '1' -> '3'
[2] '2' -> '1'
[3] '2' -> '3'
[4] '4' -> '1'


class(g1)
[1] "igraph"


png("110214.g1.png")
plot.igraph(g1)
dev.off()



bioperlを使ってて、塩基配列のアミノ酸への翻訳に成功した!キタ――o(・∀・`o)(o`・∀・´o)(o´・∀・)o キタ――♪

ちょっと時期早尚な感は否めないが、bioperlを用いてプログラミングを書いて、入力した塩基配列をアミノ酸に翻訳させてみた!!!!



$ sudo apt-get install bioperl

$ vim 110213.5.pl

a



use strict;
use warnings;
use Bio::Seq;

$naseq = Bio::Seq->new(-seq=>$ARGV[0]);
$aaseq = $naseq->translate();

print $aaseq->seq(), "\n";



Esc


$ perl 110213.5.pl "atgcatgcatgcatgc"
MHACM


キタ――o(・∀・`o)(o`・∀・´o)(o´・∀・)o キタ――♪

キタキタ━(ノ゚∀゚)ノ ┫:。・:*:・゚'★,。・:*:♪・゚'☆━━━!!!

Perlで条件文

$ vim 110213.pl

a

use strict;

use warnings;


print "好きな遺伝子を選んでください。\n";

print "1 = capase\n";

print "2 = hox\n";

print "3 = rhodopsin\n";

my $line = <STDIN>;                     # ユーザから1行入力

if ($line == 1) {                       # 1番か?

    print "アポトーシスの実行に関与します。\n";

} elsif ($line == 2) {                  # 2番か?

    print "発生時の体節の形勢に関与します。\n";

} elsif ($line == 3) {                  # 3番か?

    print "視細胞の働きを助けます\n";

} else {                                # それ以外か?

    print "ちゃんと言われた遺伝子を選んでくださいまし。\n";

}

Esc

:wq

$ perl 110213.pl
好きな遺伝子を選んでください。
1 = capase
2 = hox
3 = rhodopsin
1
アポトーシスの実行に関与します。
$ perl 110213.pl
好きな遺伝子を選んでください。
1 = capase
2 = hox
3 = rhodopsin
2
発生時の体節の形勢に関与します。
$ perl 110213.pl
好きな遺伝子を選んでください。
1 = capase
2 = hox
3 = rhodopsin
3
視細胞の働きを助けます
$ perl 110213.pl
好きな遺伝子を選んでください。
1 = capase
2 = hox
3 = rhodopsin
4
ちゃんと言われた遺伝子を選んでくださいまし。




$ vim 110213.2.pl

a


use strict;
use warnings;

print "ヒトの染色体の本数を述べなさい。\n";
my $line = <STDIN>;                     # ユーザから1行入力
if ($line == "46") {                  
    print "その通り!!\n";   
} else {
    print "違います。\n";     
}


Esc

:wq



$ perl 110213.2.pl
ヒトの染色体の本数を述べなさい。
40
違います。
$ perl 110213.2.pl
ヒトの染色体の本数を述べなさい。
46
その通り!!
kappa@kappa-desktop:~/2011/perl$





$ vim 110213.3.pl

a

print "あなたの身長を入力してください。\n";
my $x = 171;
my $line = <STDIN>;

if ($x > $line) {
print "あなたの身長$lineは日本人男性の平均身長よりも小さいです\n";
} else {
print "あなたの身長$lineは日本人男性の平均身長よりも大きいです\n";
}


Esc

:wq


$ perl 110213.3.pl
あなたの身長を入力してください。
168
あなたの身長168
は日本人男性の平均身長よりも小さいです

$ perl 110213.3.pl
あなたの身長を入力してください。
172
あなたの身長172
は日本人男性の平均身長よりも大きいです




$ vim 110214.4.pl

a

use strict;
use warnings;


print "ただいまの時刻を入力してちょ。\n";
my $line = <STDIN>;

if ($line < 8) {
print "おはよ!\n";
} elsif ($line == 12) {
print "おなかすいたっさ。\n";
} elsif ($line < 17) {
print "こんにちは!\n";
} else {
print "こんばんは。";
}

Esc


$ perl 110213.4.pl
ただいまの時刻を入力してちょ。
16
こんにちは!

Rで数式を描写してみた♪高校時代を思い出す絵になった☆

png("110212_graphics.png")

x <- seq(-5,5,length=10001)
plot(c(-5,5),c(-100,100),type="n", main = "ごく簡単な関数")
par(pch=".")


points(x,x, col = "black")
points(x,x^2, col="red")
points(x,x^3, col ="yellow")
points(x,-x,col = "green")
points(x, x^2-x , col = "blue")
points(x, x^5-x^4+x^3 , col = "grey")

dev.off()

igraph(R)によりランダムネットワークを書いてみた!!!!

$ R


install.packages("igraph")

library(igraph)

# 乱数から成る100×100のデータ行列を作成
mat <- matrix(runif(10000), nr=100)

png("110212_graph.png")
par(mfrow=c(2,1))
plot(graph.adjacency(mat>=0.9, mode="undirected"))
plot(graph.adjacency(mat>=0.9, mode="undirected"), layout=layout.circle)
dev.off()



Open Babelによる分子構造エネルギーの計算

Open bABERLという化学構造ファイル形式の変換を目的に開発されたソフトをまずはインストールする。


$ sudo apt-get install openbabel

$ obenergy CID_2244.sdf

A T O M   T Y P E S

IDX    TYPE
1    0800
2    0800
3    0801
4    0801
5    0603
6    0603
7    0603
8    0603
9    0603
10    0603
11    0601
12    0601
13    0600
14    0100
15    0100
16    0100
17    0100
18    0100
19    0100
20    0100
21    0100

C H A R G E S

IDX    CHARGE
1    0.000000
2    -0.250000
3    -0.100000
4    -0.100000
5    0.000000
6    0.000000
7    -0.100000
8    -0.100000
9    -0.100000
10    -0.100000
11    0.100000
12    0.100000
13    0.000000
14    0.100000
15    0.100000
16    0.100000
17    0.100000
18    0.000000
19    0.000000
20    0.000000
21    0.250000

S E T T I N G   U P   C A L C U L A T I O N S

SETTING UP BOND CALCULATIONS...
SETTING UP ANGLE CALCULATIONS...
SETTING UP TORSION CALCULATIONS...
SETTING UP VAN DER WAALS CALCULATIONS...
SETTING UP ELECTROSTATIC CALCULATIONS...

E N E R G Y

     TOTAL BOND STRETCHING ENERGY = 14203.861 kJ/mol
     TOTAL ANGLE BENDING ENERGY = 1685.192 kJ/mol
     TOTAL TORSIONAL ENERGY =    6.883 kJ/mol
     TOTAL VAN DER WAALS ENERGY = 17385.054 kJ/mol
     TOTAL ELECTROSTATIC ENERGY =   -3.713 kJ/mol

TOTAL ENERGY = 33277.277 kJ/mol
==============================
*** Open Babel Warning  in ReadMolecule
  WARNING: Problems reading a MDL file
Cannot read creator/dimension line line


#こういう物理化学的な計算が瞬時にできるのがLinuxを使っているメリットのひとつやろうなーー♪

mol2psによりsdf->psの変換をしたどーーーー♪



pubchemにいって、アスピリンを入力して構造式のデータをとってくる。

このとき、

1 Download SDF
2 2D SDF save

の2段階で行った。

でも、以下の用にしてもできてしまうことを発見した

$ wget http://pubchem.ncbi.nlm.nih.gov/summary/summary.cgi?cid=2244


一般化すると


$ wget http://pubchem.ncbi.nlm.nih.gov/summary/summary.cgi?cid=xxxx

で好きな化合物のsdfファイルがとってこれそうだ!!!





それでダウンロードしたファイルをカレントディレクトリに移動しといて、

$ head CID_2244.sdf
2244
  -OEChem-02121112262D

 21 21  0     0  0  0  0  0  0999 V2000
    3.7320   -0.0600    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
    6.3301    1.4400    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
    4.5981    1.4400    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
    2.8660   -1.5600    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
    4.5981   -0.5600    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.4641   -0.0600    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0


ここで、mol2psという、MOLファイル表示をPostscriptに出力するソフトウェアをインストールする。

$ wget http://merian.pch.univie.ac.at/pch/download/chemistry\
> /mol2ps/bin/mol2ps-latest-linux-i386.gz


$ gunzip mol2ps-latest-linux-i386.gz
$ sudo install mol2ps-latest-linux-i386 /usr/local/bin/mol2ps


よし、実際ps形式に変更してみるぞ!!

$ mol2ps CID_2244.sdf > CID_2244.ps


$ head CID_2244.ps
%!PS-Adobe-2.0
%%Creator: mol2ps 0.2,  Norbert Haider, University of Vienna, 2011
%%Title: CID_2244.sdf
% the following settings were used:
% font: Helvetica 14 pt (9 pt for subscripts)
% line width: 1.0
% automatic rotation:
%      0.00� around X axis
%      0.00� around Y axis
%      0.00� around Z axis

やはりブログにはpngファイルしか貼れないのでファイル形式を変換していく。


$ sudo apt-get install gv
$ ps2pdf CID_2244.ps
$ convert CID_2244.pdf -CID_2244.png






BKChemをインストール&使ってみた♪

$ sudo apt-get install bkchem

$ bkchem



局所配列・大域配列

大域的配列

2つの配列の全体を比較し、たがいに類似している領域を可能なかぎり長く、そしてギャップをできるだけ短くするように整列させるもの。Neeedleman-Wunschのアルゴリズムが基本。Needleman-Wunschのアルゴリズムは、動的計画法によってあらかじめ与えられた一致・不一致、ギャップのスコアを元に整列配列のスコアを計算し、もっとも高いスコアになるように2本の配列をペアワイズで整列させるアルゴリズム。整列のすべての可能性を探索し、その中から最適な整列を得ることができる。


$ needle

Needleman-Wunsch global alignment of two sequences
Input sequence: refseqp:NP_203124
Second sequence(s): refseqp:NP_001018443
Gap opening penalty [10.0]:    #Enter
Gap extension penalty [0.5]:   #Enter
Output alignment [np_203124.needle]:  #Enter

#生成したファイル"np_203124.needle"をのぞいてみる♪♪

$ cat np_203124.needle
########################################
# Program: needle
# Rundate: Sun 13 Feb 2011 01:05:46
# Commandline: needle
#    -asequence refseqp:NP_203124
#    -bsequence refseqp:NP_001018443
# Align_format: srspair
# Report_file: np_203124.needle
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: NP_203124
# 2: NP_001018443
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 345
# Identity:     204/345 (59.1%)
# Similarity:   241/345 (69.9%)
# Gaps:          38/345 (11.0%)
# Score: 1040.0
#
#
#=======================================

NP_203124          1 MDCVGWPPGRKWHLEKNTSCGGSSGICASYVTQMADDQGCIEEQGVEDS-     49
                                              :|......:..:...:||.  :||
NP_001018443       1 -------------------------MCRVDKEALTSENEVLEED--QDSY     23

NP_203124         50 ANEDSVDAKPDRSSFVPSLFSKKKKN---VTMRSIKTTRDRV--PTYQYN     94
                     ..||..||||||.... .||...|||   :..:....:..|:  ||:||.
NP_001018443      24 GEEDVTDAKPDRKGRF-RLFGNFKKNDGKLQEKGESESHYRIVSPTFQYK     72

NP_203124         95 MNFEKLGKCIIINNKNFDKVTGMGVRNGTDKDAEALFKCFRSLGFDVIVY    144
                     |:.:::||||||||||||:.|||.||||||:||..|||||:||||||.||
NP_001018443      73 MSHQRVGKCIIINNKNFDEKTGMNVRNGTDRDAGELFKCFKSLGFDVAVY    122

NP_203124        145 NDCSCAKMQDLLKKASEEDHTNAACFACILLSHGEENVIYGKDGVTPIKD    194
                     ||.:|..|:.|||..|||||::::||||||||||||.:|||.||..|||.
NP_001018443     123 NDQTCRNMERLLKAVSEEDHSDSSCFACILLSHGEEGMIYGTDGAMPIKT    172

NP_203124        195 LTAHFRGDRCKTLLEKPKLFFIQACRGTELDDGIQADSGPIND---TDAN    241
                     :|:.|:||.||:|:.||||||||||||:|.|||:|.||||.||   ||||
NP_001018443     173 MTSLFKGDVCKSLVGKPKLFFIQACRGSEFDDGVQTDSGPPNDTIETDAN    222

NP_203124        242 PRYKIPVEADFLFAYSTVPGYYSWRSPGRGSWFVQALCSILEEHGKDLEI    291
                     ||:||||||||||||||||||||||:||||||||||||::|.|.||.|||
NP_001018443     223 PRHKIPVEADFLFAYSTVPGYYSWRNPGRGSWFVQALCNVLSEFGKQLEI    272

NP_203124        292 MQILTRVNDRVARHFESQSDDPHFHEKKQIPCVVSMLTKELYFSQ    336
                     ||||||||..||..|||.|:||.|.||||||||||||||||||:
NP_001018443     273 MQILTRVNYMVATSFESWSEDPRFSEKKQIPCVVSMLTKELYFN-    316


#---------------------------------------
#---------------------------------------



次に、局所的整列をさせてみる。
局所的整列とは、配列の局所的に類似している領域を探し、整列させるもの、局所整列をさせるためにはSmith-Watermanのアルゴリズムが基本になる。これは上のNeedleman-Wunschのアルゴリズムの特殊な形である。Smith-Watermanのアルゴリズムでは局所的な類似性を見て整列させたいので、スコアがマイナスになるとその時点でスコアを0にしてしまう。Needleman-Wunschのアルゴリズムではは率の端から端までを整列させたが、Smith-Watermanのアルゴリズムでは、もっともスコアが高い部分から整列させる。



$ water
Smith-Waterman local alignment of sequences
Input sequence: refseqp:NP_203124
Second sequence(s): refseqp:NP_001018443
Gap opening penalty [10.0]:      #Enter
Gap extension penalty [0.5]:     #Enter
Output alignment [np_203124.water]:  #Enter


#できたファイルをのぞいてみる。

$ cat np_203124.water
########################################
# Program: needle
# Rundate: Sun 13 Feb 2011 01:05:46
# Commandline: needle
#    -asequence refseqp:NP_203124
#    -bsequence refseqp:NP_001018443
# Align_format: srspair
# Report_file: np_203124.needle
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: NP_203124
# 2: NP_001018443
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 345
# Identity:     204/345 (59.1%)
# Similarity:   241/345 (69.9%)
# Gaps:          38/345 (11.0%)
# Score: 1040.0
#
#
#=======================================

NP_203124          1 MDCVGWPPGRKWHLEKNTSCGGSSGICASYVTQMADDQGCIEEQGVEDS-     49
                                              :|......:..:...:||.  :||
NP_001018443       1 -------------------------MCRVDKEALTSENEVLEED--QDSY     23

NP_203124         50 ANEDSVDAKPDRSSFVPSLFSKKKKN---VTMRSIKTTRDRV--PTYQYN     94
                     ..||..||||||.... .||...|||   :..:....:..|:  ||:||.
NP_001018443      24 GEEDVTDAKPDRKGRF-RLFGNFKKNDGKLQEKGESESHYRIVSPTFQYK     72

NP_203124         95 MNFEKLGKCIIINNKNFDKVTGMGVRNGTDKDAEALFKCFRSLGFDVIVY    144
                     |:.:::||||||||||||:.|||.||||||:||..|||||:||||||.||
NP_001018443      73 MSHQRVGKCIIINNKNFDEKTGMNVRNGTDRDAGELFKCFKSLGFDVAVY    122

NP_203124        145 NDCSCAKMQDLLKKASEEDHTNAACFACILLSHGEENVIYGKDGVTPIKD    194
                     ||.:|..|:.|||..|||||::::||||||||||||.:|||.||..|||.
NP_001018443     123 NDQTCRNMERLLKAVSEEDHSDSSCFACILLSHGEEGMIYGTDGAMPIKT    172

NP_203124        195 LTAHFRGDRCKTLLEKPKLFFIQACRGTELDDGIQADSGPIND---TDAN    241
                     :|:.|:||.||:|:.||||||||||||:|.|||:|.||||.||   ||||
NP_001018443     173 MTSLFKGDVCKSLVGKPKLFFIQACRGSEFDDGVQTDSGPPNDTIETDAN    222

NP_203124        242 PRYKIPVEADFLFAYSTVPGYYSWRSPGRGSWFVQALCSILEEHGKDLEI    291
                     ||:||||||||||||||||||||||:||||||||||||::|.|.||.|||
NP_001018443     223 PRHKIPVEADFLFAYSTVPGYYSWRNPGRGSWFVQALCNVLSEFGKQLEI    272

NP_203124        292 MQILTRVNDRVARHFESQSDDPHFHEKKQIPCVVSMLTKELYFSQ    336
                     ||||||||..||..|||.|:||.|.||||||||||||||||||:
NP_001018443     273 MQILTRVNYMVATSFESWSEDPRFSEKKQIPCVVSMLTKELYFN-    316


#---------------------------------------
#---------------------------------------






今回のアライメントでは結果が一緒であったが、一般的に、waterでやると、無理して配列しないため、アライメントの対処から外れてしむ配列が生じがちである。
今回は、vertebrate同士の比較的近いアミノ酸配列に対してのアライメントだったからよかったけどね。

必見!BioRubyのUbuntuへのインストール方法。

RubyはデフォルトでUbuntuに入っているけど、BioRubyは入っていない。その上、apt-get installコマンドではインストールできひんかった。そこで、ソースコードからインストールを試みてみた。

$ wget http://bioruby.open-bio.org/archive/bioruby-1.4.0

#これで、http://bioruby.open-bio.org/archiveにおいてある、bioruby-1.4.0というパッケージがとってこれるぜい♪CUIでやるなら以下にプリントスクリーンを貼っておきました☆

$ cd bioruby-1.4.0

$ tar xfz bioruby-1.4.0  #これで解答できてワンサカ、カレントディレクトリにファイルが生成するはず♪


$ pwd
/home/kappa/2011/1102/bioruby-1.4.0

$ ./configure
bash: ./configure: No such file or directory

#ウギャ!なんかようわからんことになった。



#ググッた結果以下。

ーーーーーーここからーーーーーー

bash: ./configure: No such file or directory
この意味は「./」の場所に「configure」なんてファイルもディレクトリも無いですよという意味です。「./」は通常自分がいるディレクトリの場所を表しています。
今自分がいるディレクトリに(通常ソースを解凍して出来たディレクトリ)にconfigureというファイルは無いでしょうか?
無いならばhir0さんがおっしゃっているように最初からconfigureファイルがない場合もあります。

その場合は大抵ソースを解凍した中にINSTALL、README、TODOなどのファイルがあるのでまずそれを読んでみましょう。
また解凍した中にbootstrap、bootstrap.sh、autogen.sh等のファイルがありconfigure(configure.inやconfigure.acではなく純粋にただのconfigureです)が
見当たらない場合はまず最初にbootstrapやautogen.shを実行させてconfigureファイルを自動作成する必要があります。

ーーーーーここまでーーーーーー

#なるほど。全部、configure入力しまくればいいわけではないわけね。よく勉強になりました。

#README読んでちょっと頑張ってみた。

$ ruby setup.rb

$ sudo ruby setup.rb
$ ruby setup.rb config
$ ruby setup.rb install

#これでインストールできた!!・・・・はず。。。

#テスト用のコマンドもちゃんとありました♪

$ ruby setup.rb test
Running tests...
Loaded suite test
Started
....................................................................F............F...............................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................FF.................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................
Finished in 450.006201 seconds.

  1) Failure:
test_esearch_retmax_retstart(Bio::FuncTestPubmedClassMethod) [./test/functional/bio/io/test_pubmed.rb:57]:
The failure may be caused by changes of NCBI PubMed.
<35> expected but was
<20>.

  2) Failure:
test_retrieve_1id_1db(Bio::FuncTestTogoWSREST) [./test/functional/bio/io/test_togows.rb:142]:
<false> is not true.

  3) Failure:
test_libxml(Bio::TestPhyloXMLWriter_Check_LibXML) [./test/unit/bio/db/test_phyloxml_writer.rb:31]:
Error: libxml-ruby library is not present. Please install libxml-ruby library. It is needed for Bio::PhyloXML module. Unit test for PhyloXML will not be performed.
<nil> is not true.

  4) Failure:
test_libxml(Bio::TestPhyloXML_Check_LibXML) [./test/unit/bio/db/test_phyloxml.rb:29]:
Error: libxml-ruby library is not present. Please install libxml-ruby library. It is needed for Bio::PhyloXML module. Unit test for PhyloXML will not be performed.
<nil> is not true.

2564 tests, 19598 assertions, 4 failures, 0 errors




#errorは一つもなく、failureが4個出たようだ。これらに対してどう対処すればいいのか現段階ではよくわからない。


#とりあえず、気をとりなおして、BioRubyシェルを開いてみる♪♪♪

$ bioruby
Creating directory (/home/kappa/.bioruby/shell/session) ... done
Creating directory (/home/kappa/.bioruby/shell/plugin) ... done
Creating directory (/home/kappa/.bioruby/data) ... done

. . . B i o R u b y   i n   t h e   s h e l l . . .

  Version : BioRuby 1.4.0 / Ruby 1.8.7

bioruby> exit


. . . B i o R u b y   i n   t h e   s h e l l . . .

Saving history (/home/kappa/.bioruby/shell/session/history) ... done
Saving object (/home/kappa/.bioruby/shell/session/object) ... done
Saving config (/home/kappa/.bioruby/shell/session/config) ... done

#なにやら、わんさかファイルが保存されてしまっているようだ。これは解決しなければならない。biorubyを立ち上げて、閉じるつどいらないファイルが増えたんじゃかなわない笑




ちなみに、BioRubyシェルでは、BioRubyライブラリの全機能をインタラクティブに利用できるだけでなく、BioRubyの中でよく使う機能をより短いコマンドとして書くこともできるようになっている。








BLASTの置換行列を取ってきた☆

二つの配列があるとして、それらをアライメントするときは、あらかじめスコアリング方法を決めておく必要がある。これに用いられるのが置換行列である。アミノ酸や核酸の置換のしやすさを表した行列で、さまざまな置換行列が提案されている。アミノ酸の置換行列としては、BLOSUMが有名である。BLOCKSデータベース中のたんぱく質ファミリーでよく保存されている領域でアミノ酸の相対頻度と置換確率を測定し20種類のアミノ酸について全ペア間の置換それぞれについて対数オッズ比を算出したものである。BLOSUM62は、BLASTのデフォルトの置換行列であり、62%以上の配列一致性を示すアミノ酸配列間で観測された置換を基に算出して得られている。


$ wget ftp://ftp.ncbi.nih.gov/blast/matrices/BLOSUM62
$ cat BLOSUM62
#  Matrix made by matblas from blosum62.iij
#  * column uses minimum score
#  BLOSUM Clustered Scoring Matrix in 1/2 Bit Units
#  Blocks Database = /data/blocks_5.0/blocks.dat
#  Cluster Percentage: >= 62
#  Entropy =   0.6979, Expected =  -0.5209
   A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  B  Z  X  *
A  4 -1 -2 -2  0 -1 -1  0 -2 -1 -1 -1 -1 -2 -1  1  0 -3 -2  0 -2 -1  0 -4
R -1  5  0 -2 -3  1  0 -2  0 -3 -2  2 -1 -3 -2 -1 -1 -3 -2 -3 -1  0 -1 -4
N -2  0  6  1 -3  0  0  0  1 -3 -3  0 -2 -3 -2  1  0 -4 -2 -3  3  0 -1 -4
D -2 -2  1  6 -3  0  2 -1 -1 -3 -4 -1 -3 -3 -1  0 -1 -4 -3 -3  4  1 -1 -4
C  0 -3 -3 -3  9 -3 -4 -3 -3 -1 -1 -3 -1 -2 -3 -1 -1 -2 -2 -1 -3 -3 -2 -4
Q -1  1  0  0 -3  5  2 -2  0 -3 -2  1  0 -3 -1  0 -1 -2 -1 -2  0  3 -1 -4
E -1  0  0  2 -4  2  5 -2  0 -3 -3  1 -2 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4
G  0 -2  0 -1 -3 -2 -2  6 -2 -4 -4 -2 -3 -3 -2  0 -2 -2 -3 -3 -1 -2 -1 -4
H -2  0  1 -1 -3  0  0 -2  8 -3 -3 -1 -2 -1 -2 -1 -2 -2  2 -3  0  0 -1 -4
I -1 -3 -3 -3 -1 -3 -3 -4 -3  4  2 -3  1  0 -3 -2 -1 -3 -1  3 -3 -3 -1 -4
L -1 -2 -3 -4 -1 -2 -3 -4 -3  2  4 -2  2  0 -3 -2 -1 -2 -1  1 -4 -3 -1 -4
K -1  2  0 -1 -3  1  1 -2 -1 -3 -2  5 -1 -3 -1  0 -1 -3 -2 -2  0  1 -1 -4
M -1 -1 -2 -3 -1  0 -2 -3 -2  1  2 -1  5  0 -2 -1 -1 -1 -1  1 -3 -1 -1 -4
F -2 -3 -3 -3 -2 -3 -3 -3 -1  0  0 -3  0  6 -4 -2 -2  1  3 -1 -3 -3 -1 -4
P -1 -2 -2 -1 -3 -1 -1 -2 -2 -3 -3 -1 -2 -4  7 -1 -1 -4 -3 -2 -2 -1 -2 -4
S  1 -1  1  0 -1  0  0  0 -1 -2 -2  0 -1 -2 -1  4  1 -3 -2 -2  0  0  0 -4
T  0 -1  0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1  1  5 -2 -2  0 -1 -1  0 -4
W -3 -3 -4 -4 -2 -2 -3 -2 -2 -3 -2 -3 -1  1 -4 -3 -2 11  2 -3 -4 -3 -2 -4
Y -2 -2 -2 -3 -2 -1 -2 -3  2 -1 -1 -2 -1  3 -3 -2 -2  2  7 -1 -3 -2 -1 -4
V  0 -3 -3 -3 -1 -2 -2 -3 -3  3  1 -2  1 -1 -2 -2  0 -3 -1  4 -3 -2 -1 -4
B -2 -1  3  4 -3  0  1 -1  0 -3 -4  0 -3 -3 -2  0 -1 -4 -3 -3  4  1 -1 -4
Z -1  0  0  1 -3  3  4 -2  0 -3 -3  1 -1 -3 -1  0 -1 -3 -2 -2  1  4 -1 -4
X  0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2  0  0 -2 -1 -1 -1 -1 -1 -4
* -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4 -4  1








二つのアミノ酸配列に対してドットプロットを描写する。ps,pdf,tifの変換も勉強したよ♪

caspase7のヒトとzebrafishのホモログのペアワイズ作成においてドットプロットを描写する。両方の配列が同じアミノ酸や延期のときにに点が打たれる。



$ dottup
Displays a wordmatch dotplot of two sequences


Input sequence: refseqp:NP_203124    #ひとのCASP7
Second sequence: refseqp:NP_001018443   #ゼブラフィッシュのcasp7
Word size [10]:     #ここはEnter
Graph type [x11]: ps
Created dottup.ps   #psファイル作成をしている。
$ ls
110208.txt  110211.pl  110211.txt  110211_2.pl  dottup.ps  var0.pl  var1.pl


$ sudo apt-get install gv

$ gv dottup.ps
$ ps2pdf dottup.ps       #PostScriptファイルは利便性が低いのでpdfファイルに変更する。
$ ls
110208.txt  110211.pl  110211.txt  110211_2.pl  dottup.pdf  dottup.ps  var0.pl  var1.pl

$ sudo apt-get install xpdf-reader
$ xpdf dottup.pdf


$ convert dottup.pdf -dottup.tiff        #さらにこっから、pdf -> tiff 変換を行う!!!!
$ convert dottup.pdf -dottup.png      #pngファイルじゃないとブログに貼れなかったんで。。。



perlで変数2

$ vim 110211.pl

a

my $name = "programming Lesson";
print "\$name is $name\n";

my $name = "Programming" . "Lesson";
print "\$name is $name\n";

my $name = "Programming ". "Lesson";
print "\$name is $name\n";

my $name1 = "Programming";
my $name2 = "Lesson";
my $name = "$name1 $name2";
print "\$name is $name\n";

Esc

:wq


$ perl -cw 110211.pl
110211.pl syntax OK

$ perl 110211.pl #作ったプログラム110211.plの実行

$name is programming Lesson
$name is ProgrammingLesson
$name is Programming Lesson
$name is Programming Lesson

$ vim 110211_2.pl

a


use strict;
use warnings;

my $x ="123";
my $y ="456";
print $x + $y

print $x . $y;

Esc

perl -cw 110211_2.pl
1110211_2.pl syntax OK

$ perl 110211_2.pl
579

maSigProパッケージによる時系列データの処理方法

source("http://bioconductor.org/biocLite.R")
biocLite("maSigPro")

library(maSigPro)

data(edesign.abiotic)
data(data.abiotic)



edesign.abiotic
              Time Replicate Control Cold Heat Salt
Control_3H_1     3         1       1    0    0    0
Control_3H_2     3         1       1    0    0    0
Control_3H_3     3         1       1    0    0    0
Control_9H_1     9         2       1    0    0    0
Control_9H_2     9         2       1    0    0    0
Control_9H_3     9         2       1    0    0    0
Control_27H_1   27         3       1    0    0    0
Control_27H_2   27         3       1    0    0    0
Control_27H_3   27         3       1    0    0    0
Cold_3H_1        3         4       0    1    0    0
Cold_3H_2        3         4       0    1    0    0
Cold_3H_3        3         4       0    1    0    0
Cold_9H_1        9         5       0    1    0    0
Cold_9H_2        9         5       0    1    0    0
Cold_9H_3        9         5       0    1    0    0
Cold_27H_1      27         6       0    1    0    0
Cold_27H_2      27         6       0    1    0    0
Cold_27H_3      27         6       0    1    0    0
Heat_3H_1        3         7       0    0    1    0
Heat_3H_2        3         7       0    0    1    0
Heat_3H_3        3         7       0    0    1    0
Heat_9H_1        9         8       0    0    1    0
Heat_9H_2        9         8       0    0    1    0
Heat_9H_3        9         8       0    0    1    0
Heat_27H_1      27         9       0    0    1    0
Heat_27H_2      27         9       0    0    1    0
Heat_27H_3      27         9       0    0    1    0
Salt_3H_1        3        10       0    0    0    1
Salt_3H_2        3        10       0    0    0    1
Salt_3H_3        3        10       0    0    0    1
Salt_9H_1        9        11       0    0    0    1
Salt_9H_2        9        11       0    0    0    1
Salt_9H_3        9        11       0    0    0    1
Salt_27H_1      27        12       0    0    0    1
Salt_27H_2      27        12       0    0    0    1
Salt_27H_3      27        12       0    0    0    1

class(edesign.abiotic)
[1] "matrix"

data.abiotic[1,]
        Control_3H_1 Control_3H_2 Control_3H_3 Control_9H_1 Control_9H_2
STMDF90    0.1373571   -0.3653065   -0.1532945    0.4475454    0.2874768
        Control_9H_3 Control_27H_1 Control_27H_2 Control_27H_3 Cold_3H_1
STMDF90    0.2488187     0.1793259     0.1279931    -0.1173468  1.056555
        Cold_3H_2 Cold_3H_3 Cold_9H_1 Cold_9H_2 Cold_9H_3 Cold_27H_1 Cold_27H_2
STMDF90 0.3948921 0.3884203 0.8910656 0.9550419 0.8122386  0.9498117  0.7211795
        Cold_27H_3 Heat_3H_1  Heat_3H_2  Heat_3H_3 Heat_9H_1  Heat_9H_2
STMDF90  0.6432118 0.1977278 0.08401729 -0.1252850 0.3500279 0.05356246
          Heat_9H_3 Heat_27H_1 Heat_27H_2 Heat_27H_3 Salt_3H_1  Salt_3H_2
STMDF90 -0.05703404  0.1425164  0.2824874 0.03085787 0.2681767 0.06403428
        Salt_3H_3 Salt_9H_1 Salt_9H_2 Salt_9H_3 Salt_27H_1 Salt_27H_2
STMDF90 0.1757832 0.1050698 0.6643922 0.5167695  0.4754657  0.3387966
        Salt_27H_3
STMDF90  0.3596021
> class(data.abiotic)
[1] "data.frame"

#簡単のため、データをちょっといじらせてもらいます。

edesign <- edesign.abiotic[1:18, 1:4]

edesign
              Time Replicate Control Cold
Control_3H_1     3         1       1    0
Control_3H_2     3         1       1    0
Control_3H_3     3         1       1    0
Control_9H_1     9         2       1    0
Control_9H_2     9         2       1    0
Control_9H_3     9         2       1    0
Control_27H_1   27         3       1    0
Control_27H_2   27         3       1    0
Control_27H_3   27         3       1    0
Cold_3H_1        3         4       0    1
Cold_3H_2        3         4       0    1
Cold_3H_3        3         4       0    1
Cold_9H_1        9         5       0    1
Cold_9H_2        9         5       0    1
Cold_9H_3        9         5       0    1
Cold_27H_1      27         6       0    1
Cold_27H_2      27         6       0    1
Cold_27H_3      27         6       0    1


#これで、このデータの変数は、「時間」と「Treatment」の二つになった。同様にして

data <- data.abiotic[,1:18]

colname(data)
 [1] "Control_3H_1"  "Control_3H_2"  "Control_3H_3"  "Control_9H_1"
 [5] "Control_9H_2"  "Control_9H_3"  "Control_27H_1" "Control_27H_2"
 [9] "Control_27H_3" "Cold_3H_1"     "Cold_3H_2"     "Cold_3H_3"   
[13] "Cold_9H_1"     "Cold_9H_2"     "Cold_9H_3"     "Cold_27H_1"  
[17] "Cold_27H_2"    "Cold_27H_3"  

design <- make.design.matrix(edesign, degree = 2)

#"This example has three time points , so we can consider up to a quadratic regression model(degree = 2).Lager number of time points woudld potentially allow a higher polynominal degree.

class(design)
[1] "list"

design
$dis
              ColdvsControl Time TimexCold Time2 Time2xCold
Control_3H_1              0    3         0     9          0
Control_3H_2              0    3         0     9          0
Control_3H_3              0    3         0     9          0
Control_9H_1              0    9         0    81          0
Control_9H_2              0    9         0    81          0
Control_9H_3              0    9         0    81          0
Control_27H_1             0   27         0   729          0
Control_27H_2             0   27         0   729          0
Control_27H_3             0   27         0   729          0
Cold_3H_1                 1    3         3     9          9
Cold_3H_2                 1    3         3     9          9
Cold_3H_3                 1    3         3     9          9
Cold_9H_1                 1    9         9    81         81
Cold_9H_2                 1    9         9    81         81
Cold_9H_3                 1    9         9    81         81
Cold_27H_1                1   27        27   729        729
Cold_27H_2                1   27        27   729        729
Cold_27H_3                1   27        27   729        729

$groups.vector
[1] "ColdvsControl" "Control"       "ColdvsControl" "Control"     
[5] "ColdvsControl"

$edesign
              Time Replicate Control Cold
Control_3H_1     3         1       1    0
Control_3H_2     3         1       1    0
Control_3H_3     3         1       1    0
Control_9H_1     9         2       1    0
Control_9H_2     9         2       1    0
Control_9H_3     9         2       1    0
Control_27H_1   27         3       1    0
Control_27H_2   27         3       1    0
Control_27H_3   27         3       1    0
Cold_3H_1        3         4       0    1
Cold_3H_2        3         4       0    1
Cold_3H_3        3         4       0    1
Cold_9H_1        9         5       0    1
Cold_9H_2        9         5       0    1
Cold_9H_3        9         5       0    1
Cold_27H_1      27         6       0    1
Cold_27H_2      27         6       0    1
Cold_27H_3      27         6       0    1


fit <- p.vector (data, design, Q=0.05,MT.adjust = "BH",min.obs = 3)

#min.obs : genes with less than this number of true numerical values will be excluded from the analysi.Default is 3(minimun value for a quadratic fit)

tstep <- T.fit(fit, step.method = "backward", alfa = 0.05)

#次に発現変動遺伝子を得ます。

sigs <- get.siggenes(tstep, rsq = 0.6, vars ="groups")

class(sigs)
[1] "list"

names(sigs)
[1] "sig.genes" "summary"
    Control ColdvsControl
1   STMDF90       STMDF90
2   STMIA38       STMJH42
3   STMEQ29       STMDE66
4   STMEL85       STMHZ45
5   STMGU57       STMGL58
6   STMHK85       STMIF71
7   STMHJ39       STMIA38
8   STMGB57       STMEQ29
9   STMIT31       STMDW06
10  STMEY09       STMEL85
11  STMHY91       STMEG74
12  STMHS90       STMGU57
13  STMCU02       STMDV94
14  STMGB35       STMHK85
15  STMIH90       STMDV87
16  STMCF08       STMID12
17  STMDI90       STMCV66
18  STMIG08       STMGH56
19  STMCX95       STMEJ16
20  STMCO80       STMCD46
21  STMEM80       STMJE19
22  STMCF73       STMHJ39
23  STMGC06       STMJH69
24  STMEZ88       STMGB57
25  STMER65       STMIT31
26  STMHL59       STMEZ42
27  STMIK50       STMHN16
28  STMEY29       STMEY09
29  STMDT77       STMCE01
30  STMDM56       STMDG64
31  STMDM29       STMIY82
32  STMCH02       STMFB85
33  STMCS44       STMHH10
34  STMJJ85       STMGQ20
35  STMGG37       STMCY10
36  STMCH67       STMHV34
37  STMJF53       STMHY91
38  STMIQ37       STMIV44
39  STMJN55       STMDJ72
40  STMCO78       STMHS90
41  STMCM86       STMJN05
42  STMDF13       STMCH90
43  STMET96       STMIX07
44  STMIT39       STMEF65
45  STMCK87       STMCA27
46  STMIB81       STMEG36
47  STMCU87       STMGZ67
48  STMHN19       STMEB60
49  STMJI65       STMIH67
50  STMED61       STMII96
51  STMCI43       STMCU02
52  STMHU96       STMEU58
53  STMCH79       STMGM14
54  STMDU84       STMDJ03
55  STMGO62       STMGB35
56  STMIO93       STMEH38
57  STMCB61       STMGP81
58  STMEB31       STMCS39
59  STMIW42       STMCD51
60  STMGT19       STMHY77
61  STMEG09       STMGH21
62  STMCB20       STMJP11
63  STMIW62       STMIH90
64  STMHF88       STMCE67
65  STMIX47       STMIA37
66                STMCF08
67                STMHR19
68                STMHK44
69                STMDI90
70                STMIG08
71                STMEA14
72                STMCX95
73                STMCR17
74                STMCO80
75                STMDF62
76                STMDU07
77                STMDF61
78                STMIR22
79                STMIS03
80                STMEM80
81                STMCF73
82                STMJJ41
83                STMIA39
84                STMGC06
85                STMIO60
86                STMEZ88
87                STMJJ12

-----省略します------

#まずは時間に関して有意変化がありかつ、controlvscolでも有意なSTMDF90について



STMDF90 <- data[rownames(data) == "STMDF90", ]


png("STMDF90.png")
PlotGroups(STMDF90, edesign = edesign, show.fit = T, dis = design$dis, groups.vector = design$groups.vector,main = "SDMF90")
dev.off()


#see.genes(9 performs a cluster analysis to group genes by similar profiles.The resulting clusters are then plotted in two fashions:as experiment-wide expression profiles and as by-groups profies.

pdf("ColdsvsControl.pdf")
see.genes(sigs$sig.genes$ColdvsControl, main = "ColdvsControl", dis =design$dis, cluster.method="kmeans" ,cluster.data = 1, k = 9)
dev.off()

png("ColdsvsControl.png")
par(mfrow = c(9,9))
see.genes(sigs$sig.genes$ColdvsControl, main = "ColdvsControl", show.fit = T, dis =design$dis, cluster.method="kmeans" ,cluster.data = 1, k = 9)
dev.off()

png("ColdsvsControl.test.png")
s





dev.off()








マイクロアレイデータの描写のモデルを作ってみたよん♪

$ R

arrayModel <- matrix(rep(0:1,121), 11, 11)

arrayModel
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
 [1,]    0    1    0    1    0    1    0    1    0     1     0
 [2,]    1    0    1    0    1    0    1    0    1     0     1
 [3,]    0    1    0    1    0    1    0    1    0     1     0
 [4,]    1    0    1    0    1    0    1    0    1     0     1
 [5,]    0    1    0    1    0    1    0    1    0     1     0
 [6,]    1    0    1    0    1    0    1    0    1     0     1
 [7,]    0    1    0    1    0    1    0    1    0     1     0
 [8,]    1    0    1    0    1    0    1    0    1     0     1
 [9,]    0    1    0    1    0    1    0    1    0     1     0
[10,]    1    0    1    0    1    0    1    0    1     0     1
[11,]    0    1    0    1    0    1    0    1    0     1     0


png("110211_R.matrix.png")
par(mfrow = c(1,3))
image(arrayModel, col=c(1,0),main="arrayModel matrix")
image(arrayModel, col=c(0,1),main="arrayModel matrix")
image(arrayModel, col=c(2,3),main="arrayModel matrix")
dev.off()

よりマイクロアレイのシチュエーションに近付けてみるよん!!

png("110211_R.matrix.2.png")
arrayModel <- matrix(rep(0:7,121), 11, 11)
image(arrayModel, col=c(2,3,8,7,5,11,9),main="arrayModel matrix")
dev.off()


png("110211_R.matrix.3.png")
arrayModel <- matrix(rep(0:111,111*111), 7, 111)
image(arrayModel, col=c(2,3,8,7,5,11,9),main="arrayModel matrix")
dev.off()

マイクロアレイのデータは本質的には、このようなマトリックスデータである。但し、各々のスポットは連続変数であるため、今回のモデルよりも美しい絵になるわけである。