R tips：解决monocle的plot_genes_branched_heatmap函数报错

主要观点总结

在使用monocle包进行基因表达分析时遇到报错，通过查看和修改源代码解决了问题。

关键观点总结

关键观点1:

关键观点2:

关键观点3:

正文

请到「今天看啥」查看全文

fit_model_helper ( x , modelFormulaStr = trend_formula ,

expressionFamily = expressionFamily , relative_expr = relative_expr ,

disp_func = cds@dispFitInfo [[ "blind" ]] $disp_func )

if ( is . null ( model_fits ))

expression_curve as . data . frame ( matrix ( rep ( NA , nrow ( new_data )),

nrow = 1 ))

else expression_curve as . data . frame ( responseMatrix ( list ( model_fits ),

new_data , response_type = response_type ))

colnames ( expression_curve ) row . names ( new_data )

expression_curve

}, convert_to_dense = TRUE , trend_formula = trend_formula , expressionFamily = expressionFamily ,

relative_expr = relative_expr , new_data = new_data )

3 : genSmoothCurves ( new_cds [, ], cores = cores , trend_formula = trend_formula ,

relative_expr = T , new_data = rbind ( newdataA , newdataB ))

2 : plot_genes_branched_heatmap (., branch_point = 1 , num_clusters = 4 ,

show_rownames = F , return_heatmap = T )

1 : cds [ BEAM_genes , ] %>% plot_genes_branched_heatmap ( branch_point = 1 ,

num_clusters = 4 , show_rownames = F , return_heatmap = T )

梳理调用栈，并结合各节点的函数源码可以发现报错传递是：

plot_genes_branched_heatmap -> genSmoothCurves -> responseMatrix

先看一下responseMatrix的源码，如下所示，可以发现报错点的if语句那里是有问题的。因为结合最开始的报错提示表明x@family@vfamily的长度可能不为1。

monocle:::responseMatrix
function (models, newdata = NULL, response_type = "response", 
    cores = 1) 
{
    res_list  mclapply(models, function(x) {
        if (is.null(x)) {
            NA
        }
        else {
            ### 报错点 ###
            if (x@family@vfamily %in% c("negbinomial", "negbinomial.size")) {
                predict(x, newdata = newdata, type = response_type)
            }
            else if (x@family@vfamily %in% c("uninormal")) {
                predict(x, newdata = newdata, type = response_type)
            }
            else {
                10^predict(x, newdata = newdata, type = response_type)
            }
        

请到「今天看啥」查看全文