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()
