单细胞-拟时序分析-R包SCORPIUS

 注：代码转自生信技能树，仅做学习笔记。https://mp.weixin.qq.com/s/rho-EPb804bSd8w334BNQw

准备工作， 安装及加载R包，载入数据，查看数据类型。

install.packages("SCORPIUS")
library(SCORPIUS)
data(ginhoux)
ginhoux
class(ginhoux)
str(ginhoux)
c<-ginhoux[[1]]
g<-ginhoux[[2]]

c<-ginhoux[[1]]

g<-ginhoux[[2]]

expression <- ginhoux$expression
group_name <- ginhoux$sample_info$group_name
table(group_name)
dim(expression)
expression[1:4,1:4]

该数据 245个细胞，2000个基因，然后细胞分成了3组

space <- reduce_dimensionality(expression, "spearman")
draw_trajectory_plot(space, group_name, contour = TRUE)
space <- reduce_dimensionality(expression, "spearman")
draw_trajectory_plot(space, group_name, contour = TRUE)
traj <- infer_trajectory(space)
draw_trajectory_plot(space, group_name, traj$path, contour = TRUE)

# warning: setting num_permutations to 10 requires a long time (~30min) to run!
# set it to 0 and define a manual cutoff for the genes (e.g. top 200) for a much shorter execution time.
gimp <- gene_importances(
  expression, 
  traj$time, 
  num_permutations = 10, 
  num_threads = 8, 
  ntree = 10000,
  ntree_perm = 1000
) 
gimp$qvalue <- p.adjust(gimp$pvalue, "BH", length(gimp$pvalue))
gene_sel <- gimp$gene[gimp$qvalue < .05]
expr_sel <- scale_quantile(expression[,gene_sel])

# Draw a time series heatmap
time <- traj$time
draw_trajectory_heatmap(expr_sel, time)

## Also show the progression groupings
draw_trajectory_heatmap(expr_sel, time, 
                        progression_group=group_name)

## 换一个配色
table(group_name)
draw_trajectory_heatmap(
  expr_sel, time, progression_group=group_name,
  progression_group_palette = setNames(RColorBrewer::brewer.pal(3, "Set2"),
                                       levels(group_name))
)

## 基因也可以分组
modules <- extract_modules(scale_quantile(expr_sel))

draw_trajectory_heatmap(expr_sel, time,
                        progression_group=group_name, modules=modules)