circacompare_mixed()运行报错：Cholesky分解失败与收敛问题求助

使用circacompare包的circacompare_mixed()函数批量分析5000+基因的节律特征（对照组CON vs 实验组CC，每个时间点有重复样本），通过for循环处理时出现收敛失败和奇异矩阵错误：

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 1
Warning: 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.

数据结构 #

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

运行代码 #

{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)
}

系统环境 #

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 

怀疑问题出在control=list()的参数设置，寻求合适的参数配置及解决方法。

1. 简化随机效应参数 #

当前random_params=c("phi1","alpha1","k1")同时纳入相位、振幅、周期的随机效应，而样本仅3个重复id，容易导致奇异矩阵。建议逐步减少随机效应数量：

# 仅保留相位的随机效应
control=list(grouped_params=c("phi","alpha","k"), random_params=c("phi1"))
# 或仅保留振幅的随机效应
control=list(grouped_params=c("phi","alpha","k"), random_params=c("alpha1"))

若仍报错，可尝试移除所有随机效应，先用固定效应模型验证数据节律特征：

control=list(grouped_params=c("phi","alpha","k"), random_params=NULL)

2. 提升迭代次数解决收敛失败 #

根据错误提示，在control中加入nlme控制参数，提高迭代上限：

control=list(
  grouped_params=c("phi","alpha","k"), 
  random_params=c("phi1"),
  lme_control=nlme::lmeControl(msMaxIter=1000, maxIter=500)
)

3. 数据预处理优化 #

• 对基因表达量做标准化，减少数值差异对模型收敛的影响：

t[,4:5] <- scale(t[,4:5])

• 确认时间点覆盖完整周期，且每个时间点的重复样本数量均衡。

4. 批量处理容错优化 #

加入错误捕获机制，避免单个基因报错中断整个循环：

for(i in 4:5){
  gm <- colnames(t[i])
  tryCatch({
    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")),
      period = 24,alpha_threshold = 0.05
    )
    # 后续提取参数的代码
  }, error=function(e){
    message(paste("Error processing", gm, ":", e$message))
    # 可添加缺失值填充逻辑，避免结果数据框断裂
  })
}

此外，可先用circacompare()单样本模型筛选出有节律的基因，再进行混合模型分析，减少无效计算。

5. 验证模型复杂度合理性 #

仅3个生物学重复的样本量下，过度复杂的混合模型易出现奇异矩阵问题。建议先用非混合模型circacompare()验证两组节律差异，再根据结果决定是否需要混合效应模型。

内容的提问来源于stack exchange，提问作者Kylie