Entering edit mode
Hi again,
I apologise for replying to my own post, but it helps keep track if
someone will be interested.
I analysed my data with single channel analysis in limma, according to
Chapter 9. of limma usersguide.
I changed my targets file (to make it more condensed) and removed
suffix
which
identified biological replication. So my targets looks like:
>nt_trg
SlideNumber FileName Cy3 Cy5
1 93 c_093_DH_K_vs_DH_CHex.gpr hk hc
2 104 c_104_DH_CH_vs_DH_Kex.gpr hc hk
3 116 c_116_DHK_vs_DHCHex.gpr hk hc
4 16 c_016_DH_C_vs_DH_Kex.gpr hc hk
5 94 c_094_DH_K_vs_DL_Kex.gpr hk lk
6 105 c_105_DL_K_vs_DH_Kex.gpr lk hk
7 117 c_117_DHK_vs_DLKex.gpr hk lk
8 139 c_139_DL_K_vs_DH_Kex.gpr lk hk
9 92 c_092_DL_CH_vs_DL_Kex.gpr lc lk
10 106 c_106_DL_K_vs_DL_CHex.gpr lk lc
11 118 c_118_DLCH_vs_DLKex.gpr lc lk
12 23 c_023_DL_K_vs_DL_Cex.gpr lk lc
13 95 c_095_DL_CH_vs_DH_CHex.gpr lc hc
14 107 c_107_DH_CH_vs_DL_CHex.gpr hc lc
15 119 c_119_DLCH_vs_DHCHex.gpr lc hc
16 136 c_136_DH_C_vs_DL_Cex.gpr hc lc
17 101 c_101_DL_K_vs_DH_CHex.gpr lk hc
18 103 c_103_DH_CH_vs_DL_Kex.gpr hc lk
19 121 c_121_DLK_vs_DHCHex.gpr lk hc
20 15 c_015_DH_C_vs_DL_Kex.gpr hc lk
21 100 c_100_DH_K_vs_DL_CHex.gpr hk lc
22 102 c_102_DL_CH_vs_DH_Kex.gpr lc hk
23 120 c_120_DHK_vs_DLCHex.gpr hk lc
24 140 c_140_DL_C_vs_DH_Kex.gpr lc hk
I transform it to apropriate form:
>tgr_sc=targetsA2C(nt_trg)
>tgr_sc
channel.col SlideNumber FileName Target
1.1 1 93 c_093_DH_K_vs_DH_CHex.gpr hk
1.2 2 93 c_093_DH_K_vs_DH_CHex.gpr hc
2.1 1 104 c_104_DH_CH_vs_DH_Kex.gpr hc
2.2 2 104 c_104_DH_CH_vs_DH_Kex.gpr hk
3.1 1 116 c_116_DHK_vs_DHCHex.gpr hk
3.2 2 116 c_116_DHK_vs_DHCHex.gpr hc
4.1 1 16 c_016_DH_C_vs_DH_Kex.gpr hc
4.2 2 16 c_016_DH_C_vs_DH_Kex.gpr hk
5.1 1 94 c_094_DH_K_vs_DL_Kex.gpr hk
5.2 2 94 c_094_DH_K_vs_DL_Kex.gpr lk
6.1 1 105 c_105_DL_K_vs_DH_Kex.gpr lk
6.2 2 105 c_105_DL_K_vs_DH_Kex.gpr hk
7.1 1 117 c_117_DHK_vs_DLKex.gpr hk
7.2 2 117 c_117_DHK_vs_DLKex.gpr lk
8.1 1 139 c_139_DL_K_vs_DH_Kex.gpr lk
8.2 2 139 c_139_DL_K_vs_DH_Kex.gpr hk
9.1 1 92 c_092_DL_CH_vs_DL_Kex.gpr lc
9.2 2 92 c_092_DL_CH_vs_DL_Kex.gpr lk
10.1 1 106 c_106_DL_K_vs_DL_CHex.gpr lk
10.2 2 106 c_106_DL_K_vs_DL_CHex.gpr lc
11.1 1 118 c_118_DLCH_vs_DLKex.gpr lc
11.2 2 118 c_118_DLCH_vs_DLKex.gpr lk
12.1 1 23 c_023_DL_K_vs_DL_Cex.gpr lk
12.2 2 23 c_023_DL_K_vs_DL_Cex.gpr lc
13.1 1 95 c_095_DL_CH_vs_DH_CHex.gpr lc
13.2 2 95 c_095_DL_CH_vs_DH_CHex.gpr hc
14.1 1 107 c_107_DH_CH_vs_DL_CHex.gpr hc
14.2 2 107 c_107_DH_CH_vs_DL_CHex.gpr lc
15.1 1 119 c_119_DLCH_vs_DHCHex.gpr lc
15.2 2 119 c_119_DLCH_vs_DHCHex.gpr hc
16.1 1 136 c_136_DH_C_vs_DL_Cex.gpr hc
16.2 2 136 c_136_DH_C_vs_DL_Cex.gpr lc
17.1 1 101 c_101_DL_K_vs_DH_CHex.gpr lk
17.2 2 101 c_101_DL_K_vs_DH_CHex.gpr hc
18.1 1 103 c_103_DH_CH_vs_DL_Kex.gpr hc
18.2 2 103 c_103_DH_CH_vs_DL_Kex.gpr lk
19.1 1 121 c_121_DLK_vs_DHCHex.gpr lk
19.2 2 121 c_121_DLK_vs_DHCHex.gpr hc
20.1 1 15 c_015_DH_C_vs_DL_Kex.gpr hc
20.2 2 15 c_015_DH_C_vs_DL_Kex.gpr lk
21.1 1 100 c_100_DH_K_vs_DL_CHex.gpr hk
21.2 2 100 c_100_DH_K_vs_DL_CHex.gpr lc
22.1 1 102 c_102_DL_CH_vs_DH_Kex.gpr lc
22.2 2 102 c_102_DL_CH_vs_DH_Kex.gpr hk
23.1 1 120 c_120_DHK_vs_DLCHex.gpr hk
23.2 2 120 c_120_DHK_vs_DLCHex.gpr lc
24.1 1 140 c_140_DL_C_vs_DH_Kex.gpr lc
24.2 2 140 c_140_DL_C_vs_DH_Kex.gpr hk
Next, I made design matrix
>u=unique(tgr_sc$Target)
>f=factor(tgr_sc$Target,levels=u)
>design=model.matrix(~0+f)
>colnames(design)=u
>design
hk hc lk lc
1 1 0 0 0
2 0 1 0 0
3 0 1 0 0
4 1 0 0 0
5 1 0 0 0
6 0 1 0 0
7 0 1 0 0
8 1 0 0 0
9 1 0 0 0
10 0 0 1 0
11 0 0 1 0
12 1 0 0 0
13 1 0 0 0
14 0 0 1 0
15 0 0 1 0
16 1 0 0 0
17 0 0 0 1
18 0 0 1 0
19 0 0 1 0
20 0 0 0 1
21 0 0 0 1
22 0 0 1 0
23 0 0 1 0
24 0 0 0 1
25 0 0 0 1
26 0 1 0 0
27 0 1 0 0
28 0 0 0 1
29 0 0 0 1
30 0 1 0 0
31 0 1 0 0
32 0 0 0 1
33 0 0 1 0
34 0 1 0 0
35 0 1 0 0
36 0 0 1 0
37 0 0 1 0
38 0 1 0 0
39 0 1 0 0
40 0 0 1 0
41 1 0 0 0
42 0 0 0 1
43 0 0 0 1
44 1 0 0 0
45 1 0 0 0
46 0 0 0 1
47 0 0 0 1
48 1 0 0 0
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$f
[1] "contr.treatment"
*Is it correct form my design? I see, that it simply identifies what
RNA
was hybridized on each array.
>corfit=intraspotCorrelation(nt_img_lA,design)
> corfit$consensus
[1] 0.7341876
>fit=lmscFit(nt_img_lAq,design,correlation=corfit$consensus)
I want to get contrasts "hc - hk", "lc - lk", "hc - lc", "hk - lk"
and also test effect of line and temperature. To do that I write this
command:
>contr.matrix=makeContrasts(hc-hk,lc-lk,hc-lc,hk-lk,linia=(hc+hk-lc-
lk)/2,temp=(hc+lc-hk-lk)/2,inter=(hc-lc)-(hk-lk),levels=design)
* I'm not 100% sure that it's correct.
>contr.fit=contrasts.fit(fit,contr.matrix)
>contr.fit=eBayes(contr.fit)
>wynik=decideTests(contr.fit,method="global",adjust.method="BH",p.valu
e=0.05)
>summary(wynik)
hc - hk lc - lk hc - lc hk - lk linia temp inter
-1 5865 5039 3014 2685 3931 7382 1113
0 30922 33433 37177 38480 35896 28364 40776
1 6594 4909 3190 2216 3554 7635 1492