I have encountered a problem. Because my data has multiple duplicate samples at each time point, I want to use the circacompare_mixed() in the circacompare package to compared the rhythmic characteristics of gene expression between two groups (control group and experimental group), and the number of genes reached over 5000 (only two of them were extracted here). Therefore, I used the for loop again, but the results were not obtained and there were errors and warnings:
Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'!Error in nlme.formula(model = form_group, random = randomeffects_formula, :Unable to form Cholesky decomposition: the leading minor of order 1 is not pos.def.
My data: t is my data frame.
dput(t)
structure(list(id = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3,
3, 3, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3), group = c("CON",
"CON", "CON", "CON", "CON", "CON", "CON", "CON", "CON", "CON",
"CON", "CON", "CON", "CON", "CON", "CC", "CC", "CC", "CC", "CC",
"CC", "CC", "CC", "CC", "CC", "CC", "CC", "CC", "CC", "CC"),
time = c(0, 6, 12, 18, 24, 0, 6, 12, 18, 24, 0, 6, 12, 18,
24, 0, 6, 12, 18, 24, 0, 6, 12, 18, 24, 0, 6, 12, 18, 24),
Gnai3 = c(17.51364, 17.0441, 20.73355, 16.02792, 9.2444,
20.62992, 17.28127, 29.32941, 14.5062, 16.66473, 8.196866,
8.362214, 39.94942, 16.29845, 19.80216, 13.74092, 14.08834,
14.04772, 26.00864, 17.1612, 14.19446, 12.62367, 13.53981,
43.95303, 18.80314, 15.35234, 16.63615, 23.49265, 39.86695,
16.80003), Scmh1 = c(1.151244, 1.007064, 0.8978803, 0.8906483,
1.073716, 1.133559, 0.9923863, 0.3706273, 0.9956849, 1.540639,
0.9482186, 1.105261, 0.1125411, 1.31828, 1.304933, 1.016301,
1.008043, 1.098772, 0.3921483, 1.404781, 1.214758, 1.258968,
0.8213903, 0.2767058, 1.430973, 0.9601138, 1.321019, 0.6648583,
0.3420466, 1.229267)), class = "data.frame", row.names = c(NA,
30L))
My code:
{ampm<-data.frame()
amp_1m<-data.frame()
mesor_1m<-data.frame()
phase_1m<-data.frame()
amp_2m<-data.frame()
mesor_2m<-data.frame()
phase_2m<-data.frame()
amp_pm<-data.frame()
mesor_pm<-data.frame()
phase_pm<-data.frame()
rhy_p1m<-data.frame()
rhy_p2m<-data.frame()
}
for(i in 4:5){colnames(t[i])->gm
out_im<-circacompare_mixed(x=t,col_time="time", col_outcome=gm,col_group ="group",col_id="id",control=list(grouped_params=c("phi","alpha","k"),random_params=c("phi1","alpha1","k1")),period = 24,alpha_threshold = 0.05)
out_im[[2]]$gene_symbol<-gm
parameter_im<-as.data.frame(out_im[[2]])
#rhyp1m<-slice(parameter_i,1)
#rhyp2m<-slice(parameter_i,2)
mesor1m<-slice(parameter_im,3)
mesor2m<-slice(parameter_im,4)
mesorpm<-slice(parameter_im,6)
amp1m<-slice(parameter_im,7)
amp2m<-slice(parameter_im,8)
amppm<-slice(parameter_im,10)
phase1m<-slice(parameter_im,11)
phase2m<-slice(parameter_im,12)
phasepm<-slice(parameter_im,14)
amp_1m<-rbind(amp_1m,amp1m)
mesor_1m<-rbind(mesor_1m,mesor1m)
phase_1m<-rbind(phase_1m,phase1m)
amp_2m<-rbind(amp_2m,amp2m)
mesor_2m<-rbind(mesor_2m,mesor2m)
phase_2m<-rbind(phase_2m,phase2m)
#rhy_p1m<-rbind(rhy_p1m,rhyp1m)
#rhy_p2m<-rbind(rhy_p2m,rhyp2m)
amp_pm<-rbind(amp_pm,amppm)
mesor_pm<-rbind(mesor_pm,mesorpm)
phase_pm<-rbind(phase_pm,phasepm)
}
My running results:
Warning: Iteration 1, LME step: nlminb() did not converge (code = 1). Do increase 'msMaxIter'! PORT message: false convergence (8)Warning: Singular precision matrix in level -1, block 1Warning: Singular precision matrix in level -1, block . Error in nlme.formula(model = form_group, random = randomeffects_formula, :
Unable to form Cholesky decomposition: the leading minor of order 1 is not pos.def.
My system environment:
R version 4.2.1 (2022-06-23 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)
Matrix products: default
locale:
[1] LC_COLLATE=Chinese (Simplified)_China.utf8
[2] LC_CTYPE=Chinese (Simplified)_China.utf8
[3] LC_MONETARY=Chinese (Simplified)_China.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=Chinese (Simplified)_China.utf8
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
loaded via a namespace (and not attached):
[1] fgsea_1.24.0 colorspace_2.0-3
[3] ggtree_3.6.2 rjson_0.2.21
[5] gson_0.0.9 circlize_0.4.15
[7] qvalue_2.30.0 XVector_0.38.0
[9] GlobalOptions_0.1.2 aplot_0.1.9
[11] clue_0.3-64 rstudioapi_0.14
[13] farver_2.1.1 graphlayouts_0.8.4
[15] ggrepel_0.9.2 bit64_4.0.5
[17] AnnotationDbi_1.60.0 fansi_1.0.3
[19] scatterpie_0.1.8 codetools_0.2-19
[21] splines_4.2.1 doParallel_1.0.17
[23] cachem_1.0.6 GOSemSim_2.24.0
[25] knitr_1.42 polyclip_1.10-4
[27] jsonlite_1.8.4 cluster_2.1.4
[29] GO.db_3.16.0 png_0.1-8
[31] ggforce_0.4.1 compiler_4.2.1
[33] httr_1.4.5 assertthat_0.2.1
[35] Matrix_1.5-3 fastmap_1.1.0
[37] lazyeval_0.2.2 cli_3.6.0
[39] tweenr_2.0.2 htmltools_0.5.4
[41] tools_4.2.1 igraph_1.3.5
[43] gtable_0.3.1 glue_1.6.2
[45] GenomeInfoDbData_1.2.9 reshape2_1.4.4
[47] dplyr_1.0.10 fastmatch_1.1-3
[49] Rcpp_1.0.9 enrichplot_1.18.3
[51] Biobase_2.58.0 vctrs_0.5.1
[53] Biostrings_2.66.0 ape_5.6-2
[55] nlme_3.1-161 iterators_1.0.14
[57] ggraph_2.1.0 xfun_0.36
[59] stringr_1.5.0 lifecycle_1.0.3
[61] clusterProfiler_4.6.0 DOSE_3.24.2
[63] zlibbioc_1.44.0 MASS_7.3-58.1
[65] scales_1.2.1 tidygraph_1.2.2
[67] parallel_4.2.1 RColorBrewer_1.1-3
[69] yaml_2.3.7 ComplexHeatmap_2.14.0
[71] memoise_2.0.1 gridExtra_2.3
[73] ggplot2_3.4.1 downloader_0.4
[75] ggfun_0.0.9 HDO.db_0.99.1
[77] yulab.utils_0.0.6 stringi_1.7.12
[79] RSQLite_2.2.20 circacompare_0.1.1
[81] S4Vectors_0.36.1 foreach_1.5.2
[83] tidytree_0.4.2 BiocGenerics_0.44.0
[85] BiocParallel_1.32.5 shape_1.4.6
[87] GenomeInfoDb_1.34.9 rlang_1.0.6
[89] pkgconfig_2.0.3 matrixStats_0.63.0
[91] bitops_1.0-7 evaluate_0.20
[93] lattice_0.20-45 purrr_1.0.1
[95] treeio_1.22.0 patchwork_1.1.2
[97] cowplot_1.1.1 shadowtext_0.1.2
[99] bit_4.0.5 tidyselect_1.2.0
[101] plyr_1.8.8 magrittr_2.0.3
[103] R6_2.5.1 IRanges_2.32.0
[105] generics_0.1.3 DBI_1.1.3
[107] pillar_1.8.1 withr_2.5.0
[109] KEGGREST_1.38.0 RCurl_1.98-1.9
[111] tibble_3.1.8 crayon_1.5.2
[113] utf8_1.2.2 rmarkdown_2.20
[115] viridis_0.6.2 GetoptLong_1.0.5
[117] grid_4.2.1 data.table_1.14.6
[119] blob_1.2.3 digest_0.6.31
[121] tidyr_1.2.1 gridGraphics_0.5-1
[123] stats4_4.2.1 munsell_0.5.0
[125] viridisLite_0.4.1 ggplotify_0.1.0
I guess it may be a problem with the parameters in control=list(), but I'm not sure how to choose the appropriate parameters for my data.
