monocle2也可以用于bulk RNA-seq分析吗？

rm(list = ls())  
# 这个步骤超级耗费时间，拟时序分析在单细胞领域火爆起来的
# 这个项目就一百多个普通转录组样品，跟单细胞动辄成千上万个细胞还是有差距
# 所以个人电脑仍然是可以hold住 
library(monocle)
library(plyr)
packageVersion("igraph")

early_state <- function(cds){
if (length(unique(pData(cds)$State)) > 1){
    progenitorState <- table(pData(cds)$State, pData(cds)$Grade)[,"Normal"]
    return(as.numeric(names(progenitorState)[which (progenitorState == max(progenitorState))]))
  } else {
    return (1)
  }
}

label <- read.delim('../data/monocle_label.txt', sep = '\t', header = TRUE, row.names = 1, check.names = FALSE)
head(label)

expr <- read.delim('../data/monocle_expr.txt', sep = '\t', header = TRUE, row.names = 2, check.names = FALSE)
expr$ID <- NULL
identical(sort(rownames(label)), sort(colnames(expr)))

fd <- data.frame(row.names=rownames(expr), gene_short_name=rownames(expr)); head(fd)
monocleObj <- newCellDataSet(as.matrix(expr), 
                             phenoData = new("AnnotatedDataFrame", data = label), 
                             featureData=new("AnnotatedDataFrame", data = fd), 
                             expressionFamily=negbinomial.size())
# 重点就是构建对象，后续分析全部是使用这个对象
monocleObj

pData(monocleObj)$Grade <- factor(pData(monocleObj)$Grade, levels = c('Normal', 'FibCS', 'DN', 'T1', 'T2', 'T3-4', 'Mixed'))
pData(monocleObj)$site <- factor(pData(monocleObj)$site, levels = c('Normal', 'FibCS', 'DN', 'Tumor'))
pData(monocleObj)$pretx <- factor(pData(monocleObj)$pretx, levels = c('Yes', 'No', 'Not applicable'))
pData(monocleObj)$op <- factor(pData(monocleObj)$op, levels = c('Total_hepatectomy', 'Partial_hepatectomy', 'Not applicable'))
head(pData(monocleObj))

monocleObj <- estimateSizeFactors(monocleObj)
monocleObj <- estimateDispersions(monocleObj)
monocleObj <- detectGenes(monocleObj, min_expr = 1)

head(fData(monocleObj))
expressed_genes <- row.names(subset(fData(monocleObj), num_cells_expressed >= 10))
length(expressed_genes)  # 17814 genes

disp_table <- dispersionTable(monocleObj); 
rownames(disp_table) <- disp_table$gene_id

### Dispersion genes ###
monocleObj_ordering_genes <- subset(disp_table, mean_expression >= 10 & dispersion_empirical >= 2 * dispersion_fit)$gene_id; length(monocleObj_ordering_genes) # 2127
disp_table$cols <- 'no'; disp_table[as.character(monocleObj_ordering_genes), 'cols'] <- 'in use'
ggplot(disp_table, aes(log10(mean_expression+1), dispersion_empirical, col = cols)) + 
  geom_point() +
  scale_color_manual(values = c('red2', 'grey80'))

monocleObj <- setOrderingFilter(monocleObj, ordering_genes = monocleObj_ordering_genes)
monocleObj <- reduceDimension(monocleObj, method = 'DDRTree')
monocleObj <- orderCells(monocleObj)
monocleObj <- orderCells(monocleObj, root_state = early_state(monocleObj))
head(pData(monocleObj))

plot_cell_trajectory(monocleObj, color_by = 'site') + 
  scale_color_manual(values = c('Normal'='#bdc1c8', 'FibCS'='#67814f', 'DN'='#f2ae2e', 'Tumor'='#e65251'))
ggsave('../results/disp.site.orig.pdf', units = 'cm', width = 9, height = 10)