Hello,
I have looked at all of the examples and previous posts but have been unable to find my particular issue. Basically, I am receiving an error (Error in rowMeans(unwrapdups(data[, index], ndups = ndups, spacing = spacing)) : 'x' must be an array of at least two dimensions) when I try to use limmaCTData with my qPCR set. If i change to ndups=2 (with or without reorganizing the data), the error disappears but my gene names are not conserved in the heatmapSig(). My design matrix has duplicate samples in a column but I beleive the ndups is only for feature names, correct? I hope the following is clear but if you need more info or context please let me know. Thank you!
>raw <- readCtData(files = "1131361180.csv", path = path, format = "BioMark", n.features =48, n.data = 48)
>raw.filt <- raw[, c(1:21, 23:26, 28, 29, 31:42, 44:48)]
>d.norm.raw.filt <- normalizeCtData(raw.filt, deltaCt.genes = "ACTB")
>print(d.norm.raw.filt)
An object of class "qPCRset"
Size: 48 features, 44 samples
Feature types:
Feature names: PD1 BCL6 BATF ...
Feature classes:
Feature categories: OK
Sample names: 1P3SB 1CPASB 7PICNS ...
>covariates <- pData(d.norm.raw.filt)
>design <- (model.matrix(~0 +Day:Adj:Ant, covariates))[, -c(2, 16)]
>colnames(design)<- PDataColnames
>contrasts <- makeContrasts(D1CPANS-D1NNS, D1CPBNS-D1NNS, D7CPBNS-D7NNS, D1MPNS-D1NNS, D7MPNS-D7NNS, D1P3NS-D1NNS, D7P3NS-D7NNS, D1PICNS-D1NNS, D7PICNS-D7NNS, D1R8NS-D1NNS, D7R8NS-D7NNS, D1CPASB-D1NSB, D1CPBSB-D1NSB, D7CPBSB-D7NSB, D1MPSB-D1NSB, D7MPSB-D7NSB, D1P3SB-D1NSB, D7P3SB-D7NSB, D1PICSB-D1NSB, D7PICSB-D7NSB, D1R8SB-D1NSB, D7R8SB-D7NSB, levels= design)
>qDE.limma <- limmaCtData(new.d.norm.raw.filt, design = design, contrasts = contrasts, sort = TRUE, stringent = TRUE, ndups = 1, spacing = 1)
Error in rowMeans(unwrapdups(data[, index], ndups = ndups, spacing = spacing)) :
'x' must be an array of at least two dimensions
I have tried many things to fix this and can't seem to except to make ndups=2 but if I do so my gene names are not conserved in subsequent plots. Any help would be greatly appreciated.
Thank you!
-Julio
R version 3.2.1 (2015-06-18)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.2 (Yosemite)
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] parallel stats graphics grDevices utils datasets
[7] methods base
other attached packages:
[1] statmod_1.4.21 NMF_0.20.6 cluster_2.0.3
[4] rngtools_1.2.4 pkgmaker_0.22 registry_0.3
[7] HTqPCR_1.22.0 limma_3.24.14 RColorBrewer_1.1-2
[10] Biobase_2.28.0 BiocGenerics_0.14.0
loaded via a namespace (and not attached):
[1] Rcpp_0.12.0 BiocInstaller_1.18.4
[3] plyr_1.8.3 bitops_1.0-6
[5] iterators_1.0.7 tools_3.2.1
[7] zlibbioc_1.14.0 digest_0.6.8
[9] preprocessCore_1.30.0 gtable_0.1.2
[11] gridBase_0.4-7 foreach_1.4.2
[13] proto_0.3-10 stringr_1.0.0
[15] gtools_3.5.0 caTools_1.17.1
[17] stats4_3.2.1 grid_3.2.1
[19] gdata_2.17.0 ggplot2_1.0.1
[21] reshape2_1.4.1 magrittr_1.5
[23] MASS_7.3-43 gplots_2.17.0
[25] scales_0.2.5 codetools_0.2-14
[27] xtable_1.7-4 colorspace_1.2-6
[29] KernSmooth_2.23-15 stringi_0.5-5
[31] affy_1.46.1 munsell_0.4.2
[33] doParallel_1.0.8 affyio_1.36.0
Thank you, Gordon, for your input. Unfortunately when I try to do them separately I get the following error in the lmFit function:
fit <- lmFit(d.norm.raw.filt, design = design)
Error in as.vector(data) :
no method for coercing this S4 class to a vector
Some of my samples have duplicates and some do not. When I generate my "design" matrix using model.matrix it creates columns for my samples but also additional columns in which none of my samples correspond to the column (columns 2 and 16 specifically):
design <- (model.matrix(~0 +Day:Adj:Ant, covariates))
print(design)
DayDay 1:AdjCPG-ODN-A:AntNo Stim. DayDay 7:AdjCPG-ODN-A:AntNo Stim. DayDay 1:AdjCPG-ODN-B:AntNo Stim. DayDay 7:AdjCPG-ODN-B:AntNo Stim.
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0
6 0 0 0 0
7 0 0 0 0
8 1 0 0 0
9 0 0 0 0
10 0 0 0 0
11 0 0 0 0
12 0 0 0 0
13 0 0 0 0
14 0 0 1 0
15 0 0 0 0
16 0 0 0 0
17 0 0 0 0
18 0 0 0 0
19 0 0 0 0
20 0 0 0 0
21 0 0 0 0
22 0 0 0 0
23 0 0 0 0
24 0 0 0 0
25 0 0 0 0
26 0 0 0 0
27 0 0 0 0
28 0 0 0 0
29 0 0 0 0
30 0 0 0 0
31 0 0 0 0
32 0 0 0 0
33 0 0 0 1
34 0 0 0 0
35 0 0 0 0
36 0 0 0 0
37 0 0 0 0
38 0 0 0 0
39 0 0 0 0
40 0 0 0 0
41 0 0 0 0
42 0 0 0 0
43 0 0 0 0
44 0 0 0 0
DayDay 1:AdjMPLA:AntNo Stim. DayDay 7:AdjMPLA:AntNo Stim. DayDay 1:AdjNone:AntNo Stim. DayDay 7:AdjNone:AntNo Stim. DayDay 1:AdjPam3CSK4:AntNo Stim.
1 0 0 0 0 0
2 0 0 0 0 0
3 0 0 0 0 0
4 0 0 0 0 0
5 0 0 0 0 0
6 0 0 0 0 0
7 0 0 0 0 1
8 0 0 0 0 0
9 0 0 0 0 0
10 0 0 0 0 1
11 0 0 0 0 0
12 0 1 0 0 0
13 0 0 0 0 0
14 0 0 0 0 0
15 0 1 0 0 0
16 0 0 0 0 0
17 0 0 0 0 0
18 0 0 0 0 0
19 0 0 0 0 0
20 0 0 0 0 0
21 0 0 0 0 0
22 0 0 1 0 0
23 0 0 0 0 0
24 0 0 0 0 0
25 0 0 1 0 0
26 0 0 0 0 0
27 0 0 0 0 0
28 1 0 0 0 0
29 0 0 0 0 0
30 0 0 0 0 0
31 1 0 0 0 0
32 0 0 0 0 0
33 0 0 0 0 0
34 0 0 0 0 0
35 0 0 0 0 0
36 0 0 0 0 0
37 0 0 0 0 0
38 0 0 0 0 0
39 0 0 0 0 0
40 0 0 0 0 0
41 0 0 0 1 0
42 0 0 0 0 0
43 0 0 0 0 0
44 0 0 0 1 0
DayDay 7:AdjPam3CSK4:AntNo Stim. DayDay 1:AdjPoly-I:C:AntNo Stim. DayDay 7:AdjPoly-I:C:AntNo Stim. DayDay 1:AdjR848:AntNo Stim. DayDay 7:AdjR848:AntNo Stim.
1 0 0 0 0 0
2 0 0 0 0 0
3 0 0 1 0 0
4 0 0 0 0 0
5 0 0 0 0 0
6 0 0 0 0 0
7 0 0 0 0 0
8 0 0 0 0 0
9 0 0 0 0 0
10 0 0 0 0 0
11 0 0 0 0 0
12 0 0 0 0 0
13 0 0 0 0 0
14 0 0 0 0 0
15 0 0 0 0 0
16 0 0 0 0 0
17 0 0 0 0 0
18 0 0 0 0 0
19 0 1 0 0 0
20 0 0 0 0 0
21 0 0 0 0 0
22 0 0 0 0 0
23 0 0 0 0 1
24 0 0 0 0 0
25 0 0 0 0 0
26 0 0 0 0 0
27 0 0 0 0 0
28 0 0 0 0 0
29 0 0 0 0 0
30 0 0 0 0 0
31 0 0 0 0 0
32 1 0 0 0 0
33 0 0 0 0 0
34 0 0 0 0 0
35 1 0 0 0 0
36 0 0 0 0 0
37 0 0 0 0 0
38 0 0 0 0 0
39 0 0 0 0 0
40 0 0 0 0 0
41 0 0 0 0 0
42 0 0 0 1 0
43 0 0 1 0 0
44 0 0 0 0 0
DayDay 1:AdjCPG-ODN-A:AntSEB DayDay 7:AdjCPG-ODN-A:AntSEB DayDay 1:AdjCPG-ODN-B:AntSEB DayDay 7:AdjCPG-ODN-B:AntSEB DayDay 1:AdjMPLA:AntSEB
1 0 0 0 0 0
2 1 0 0 0 0
3 0 0 0 0 0
4 0 0 0 0 0
5 1 0 0 0 0
6 0 0 0 0 0
7 0 0 0 0 0
8 0 0 0 0 0
9 0 0 0 0 0
10 0 0 0 0 0
11 0 0 1 0 0
12 0 0 0 0 0
13 0 0 0 0 0
14 0 0 0 0 0
15 0 0 0 0 0
16 0 0 0 0 0
17 0 0 0 0 0
18 0 0 0 0 0
19 0 0 0 0 0
20 0 0 0 0 0
21 0 0 0 0 0
22 0 0 0 0 0
23 0 0 0 0 0
24 0 0 0 0 1
25 0 0 0 0 0
26 0 0 0 0 1
27 0 0 0 0 0
28 0 0 0 0 0
29 0 0 0 0 0
30 0 0 0 1 0
31 0 0 0 0 0
32 0 0 0 0 0
33 0 0 0 0 0
34 0 0 0 0 0
35 0 0 0 0 0
36 0 0 0 0 0
37 0 0 0 0 0
38 0 0 0 0 0
39 0 0 0 0 0
40 0 0 0 0 0
41 0 0 0 0 0
42 0 0 0 0 0
43 0 0 0 0 0
44 0 0 0 0 0
DayDay 7:AdjMPLA:AntSEB DayDay 1:AdjNone:AntSEB DayDay 7:AdjNone:AntSEB DayDay 1:AdjPam3CSK4:AntSEB DayDay 7:AdjPam3CSK4:AntSEB DayDay 1:AdjPoly-I:C:AntSEB
1 0 0 0 1 0 0
2 0 0 0 0 0 0
3 0 0 0 0 0 0
4 0 0 0 1 0 0
5 0 0 0 0 0 0
6 1 0 0 0 0 0
7 0 0 0 0 0 0
8 0 0 0 0 0 0
9 1 0 0 0 0 0
10 0 0 0 0 0 0
11 0 0 0 0 0 0
12 0 0 0 0 0 0
13 0 0 0 0 0 1
14 0 0 0 0 0 0
15 0 0 0 0 0 0
16 0 0 0 0 0 1
17 0 1 0 0 0 0
18 0 0 0 0 0 0
19 0 0 0 0 0 0
20 0 1 0 0 0 0
21 0 0 0 0 0 0
22 0 0 0 0 0 0
23 0 0 0 0 0 0
24 0 0 0 0 0 0
25 0 0 0 0 0 0
26 0 0 0 0 0 0
27 0 0 0 0 1 0
28 0 0 0 0 0 0
29 0 0 0 0 1 0
30 0 0 0 0 0 0
31 0 0 0 0 0 0
32 0 0 0 0 0 0
33 0 0 0 0 0 0
34 0 0 0 0 0 0
35 0 0 0 0 0 0
36 0 0 1 0 0 0
37 0 0 0 0 0 0
38 0 0 0 0 0 0
39 0 0 1 0 0 0
40 0 0 0 0 0 0
41 0 0 0 0 0 0
42 0 0 0 0 0 0
43 0 0 0 0 0 0
44 0 0 0 0 0 0
DayDay 7:AdjPoly-I:C:AntSEB DayDay 1:AdjR848:AntSEB DayDay 7:AdjR848:AntSEB
1 0 0 0
2 0 0 0
3 0 0 0
4 0 0 0
5 0 0 0
6 0 0 0
7 0 0 0
8 0 0 0
9 0 0 0
10 0 0 0
11 0 0 0
12 0 0 0
13 0 0 0
14 0 0 0
15 0 0 0
16 0 0 0
17 0 0 0
18 0 0 1
19 0 0 0
20 0 0 0
21 0 0 1
22 0 0 0
23 0 0 0
24 0 0 0
25 0 0 0
26 0 0 0
27 0 0 0
28 0 0 0
29 0 0 0
30 0 0 0
31 0 0 0
32 0 0 0
33 0 0 0
34 0 1 0
35 0 0 0
36 0 0 0
37 0 1 0
38 1 0 0
39 0 0 0
40 1 0 0
41 0 0 0
42 0 0 0
43 0 0 0
44 0 0 0
attr(,"assign")
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Day
[1] "contr.treatment"
attr(,"contrasts")$Adj
[1] "contr.treatment"
attr(,"contrasts")$Ant
[1] "contr.treatment"
I eliminated two columns from my matrix given that they did not correspond to samples I actually had (they had a 0 value straight down) and then renamed them.
design <- (model.matrix(~0 +Day:Adj:Ant, covariates))[, -c(2, 16)]
Here is what design looks like afterwards (although affected by my copy pasting of course)
> colnames(design)<- PDataColnames
print(design)
D1CPANS D1CPBNS D7CPBNS D1MPNS D7MPNS D1NNS D7NNS D1P3NS D7P3NS D1PICNS D7PICNS D1R8NS D7R8NS D1CPASB D1CPBSB D7CPBSB D1MPSB D7MPSB D1NSB
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
3 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
5 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
7 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
8 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
10 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
12 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
14 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
15 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
17 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
18 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
19 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
21 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
22 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
23 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
24 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
25 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
26 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
27 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
28 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
29 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
30 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
31 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
32 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
33 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
34 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
35 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
36 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
37 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
38 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
39 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
40 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
41 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
42 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
43 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
44 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
D7NSB D1P3SB D7P3SB D1PICSB D7PICSB D1R8SB D7R8SB
1 0 1 0 0 0 0 0
2 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0
4 0 1 0 0 0 0 0
5 0 0 0 0 0 0 0
6 0 0 0 0 0 0 0
7 0 0 0 0 0 0 0
8 0 0 0 0 0 0 0
9 0 0 0 0 0 0 0
10 0 0 0 0 0 0 0
11 0 0 0 0 0 0 0
12 0 0 0 0 0 0 0
13 0 0 0 1 0 0 0
14 0 0 0 0 0 0 0
15 0 0 0 0 0 0 0
16 0 0 0 1 0 0 0
17 0 0 0 0 0 0 0
18 0 0 0 0 0 0 1
19 0 0 0 0 0 0 0
20 0 0 0 0 0 0 0
21 0 0 0 0 0 0 1
22 0 0 0 0 0 0 0
23 0 0 0 0 0 0 0
24 0 0 0 0 0 0 0
25 0 0 0 0 0 0 0
26 0 0 0 0 0 0 0
27 0 0 1 0 0 0 0
28 0 0 0 0 0 0 0
29 0 0 1 0 0 0 0
30 0 0 0 0 0 0 0
31 0 0 0 0 0 0 0
32 0 0 0 0 0 0 0
33 0 0 0 0 0 0 0
34 0 0 0 0 0 1 0
35 0 0 0 0 0 0 0
36 1 0 0 0 0 0 0
37 0 0 0 0 0 1 0
38 0 0 0 0 1 0 0
39 1 0 0 0 0 0 0
40 0 0 0 0 1 0 0
41 0 0 0 0 0 0 0
42 0 0 0 0 0 0 0
43 0 0 0 0 0 0 0
44 0 0 0 0 0 0 0
Warning messages:
1: display list redraw incomplete
2: display list redraw incomplete
3: display list redraw incomplete
All columns have at least one "1" value in them now.
So if i proceed:
fit <- lmFit(d.norm.raw.filt, design = design)
Error in as.vector(data) :
no method for coercing this S4 class to a vector
Even if I don't filter out these two columns I still get the same error. Is it my matrix?
The error message has nothing to do with the design matrix. It is rather caused by the fact that the d.norm.raw.filt data object is specific to the HTqPCR package and so isn't understood directly by limma. So you need to use exprs(d.norm.raw) to turn it into a matrix.
However your design matrix does look very complicated, and most of your groups seem to have no replication. There isn't anything to be gained from a 3-way interaction with so many groups. Consider instead using the approach explained in Section 9.3 of the limma User's Guide.
Thanks, Gordon; I read Section 9.3. My goal here is to make several comparisons, particularly information I would get from a HeatmapSig and a plotCtRQ. Do you happen to know why it is that I can get a qDE value from limmaCTData when ndups=2 and not when ndups =1? I've tried to unwrap limmaCTData as you've said but I get plenty of warnings and am ulitimately not able to use heatmapSig or plotCtRQ afterwards. I am very new to R so if some of this is obvious in some way, I apologize. I do not mind using ndups=2 in limmaCTData if this would preserve the gene names so that I would not have to go back and look at what these are when I use plotCtRQ.
fit <- lmFit(exprs(d.norm.raw.filt), design = design)
> fit2 <- contrasts.fit(fit, contrasts=contrasts)
Warning messages:
1: In data.row.names(row.names, rowsi, i) :
some row.names duplicated: 4,5,6,10,11,12,16,17,18,22,23,24,28,29,30,34,35,36,40,41,42,46,47,48 --> row.names NOT used
2: In data.row.names(row.names, rowsi, i) :
some row.names duplicated: 4,5,6,10,11,12,16,17,18,22,23,24,28,29,30,34,35,36,40,41,42,46,47,48 --> row.names NOT used
fit2<- eBayes(fit2)
Warning messages:
1: In data.row.names(row.names, rowsi, i) :
some row.names duplicated: 4,5,6,10,11,12,16,17,18,22,23,24,28,29,30,34,35,36,40,41,42,46,47,48 --> row.names NOT used
2: In data.row.names(row.names, rowsi, i) :
some row.names duplicated: 4,5,6,10,11,12,16,17,18,22,23,24,28,29,30,34,35,36,40,41,42,46,47,48 --> row.names NOT used
heatmapSig(fit2)
Error: Two subscripts required
Thanks again for your help, Gordon.
Have you tried using topTable(), as I suggested? That would immediately give you results for any of your comparisons.
The warning messages from contrast.fit() and eBayes() don't matter, but you could avoid them by averaging the replicate probes using avereps() before you run lmFit().
As far as I can tell, limmaCtData() doesn't work because it contains a bug. You can't set ndups=2 because it would be a misrepresentation of your data. Analysing PCR data is not very difficult using methods designed for microarrays so you don't strictly need the HTqPCR package. But if you insist on using the HTqPCR package (using the functions HeamapSig and plotCtRQ for example) then you will have to write to the HTqPCR authors directly and persuade them to fix the problem. No one else will be able to help you.
Note that I am not an HTqPCR author and I cannot answer questions about it.
Thank you, Gordon for all of your help! Yes, I tried the topTable function and was able to get results. Would there be a way to visualize these results? What other packages could I use to analyze my results besides the HTqPCR package? Do you know how I could obtain the email information of the HTqPCR authors? Thanks again!
Best,
-Julio
Well, there are any number of ways to plot the results. Have you tried plotMD(fit2) or plotMDS(d.norm.raw.filt)? For heatmaps, there is heatmap.2 in the gplots package.