学习小组D8

单样本和多样本拟时序～

其实并没有太多的不同，因为都是针对seurat对象的，多样本也整合过

如果要说有区别，可能多样本可以做基于treat/ctrl的比较（之后探索

谢谢花花！

#单样本拟时序
rm(list = ls())
library(Seurat)
library(monocle)
library(dplyr)
load("sce.Rdata")
scRNA = sce
table(Idents(scRNA))
head(scRNA@meta.data)
DimPlot(scRNA,label = T)
scRNA$celltype = Idents(scRNA) #celltype是细胞类型注释
levels(Idents(scRNA)) #打出来细胞类型供复制
scRNA <- subset(scRNA, idents = c(" NK", " CD8 T"))
table(Idents(scRNA))
#随机种子抽样
set.seed(1234)
scRNA = subset(scRNA,downsample = 100)

#创建monocle可以接受的CellDataSet对象
# count矩阵，官方建议用count
ct <- scRNA@assays$RNA$counts
# 基因注释
gene_ann <- data.frame(
  gene_short_name = row.names(ct), 
  row.names = row.names(ct)
)
fd <- new("AnnotatedDataFrame",
          data=gene_ann)
# 临床信息
pd <- new("AnnotatedDataFrame",
          data=scRNA@meta.data)
#新建CellDataSet对象
sc_cds <- newCellDataSet(
  ct, 
  phenoData = pd,
  featureData =fd,
  expressionFamily = negbinomial.size(),
  lowerDetectionLimit=1)
sc_cds

#构建细胞发育轨迹
sc_cds <- estimateSizeFactors(sc_cds) #估计大小因子
sc_cds <- estimateDispersions(sc_cds) #估计离散度
fdif = "diff_test_res.Rdata"
if(!file.exists(fdif)){
  diff_test_res <- differentialGeneTest(sc_cds,
                                        fullModelFormulaStr = "~celltype",
                                        cores = 4)
  save(diff_test_res,file = fdif)
}
load(fdif)
ordering_genes <- row.names(subset(diff_test_res, qval < 0.01))
#查看基因，筛选适合用于排序的，设置为排序要使用的基因
head(ordering_genes)
length(ordering_genes)
sc_cds <- setOrderingFilter(sc_cds, ordering_genes)
#画出选择的基因
plot_ordering_genes(sc_cds)
#降维
sc_cds <- reduceDimension(sc_cds)
#细胞排序
sc_cds <- orderCells(sc_cds)

#绘图展示
##发育轨迹图
library(ggsci)
p1 = plot_cell_trajectory(sc_cds)+ scale_color_nejm()
p2 = plot_cell_trajectory(sc_cds, color_by = 'Pseudotime') 
p3 = plot_cell_trajectory(sc_cds, color_by = 'celltype')  + scale_color_npg()
library(patchwork)
p2+p1/p3

##经典的拟时序热图
gene_to_cluster = diff_test_res %>% arrange(qval) %>% head(50) %>% pull(gene_short_name);head(gene_to_cluster)
plot_pseudotime_heatmap(sc_cds[gene_to_cluster,],
                        num_clusters = nlevels(Idents(scRNA)), 
                        show_rownames = TRUE,
                        cores = 4,return_heatmap = TRUE,
                        hmcols = colorRampPalette(c("navy", "white", "firebrick3"))(100))

##num_clusters = nlevels(Idents(scRNA))是把基因聚成几簇，推荐是由几类细胞就设置几，用nlevels实现了自动化（花花原创）。
#colorRampPalette(c("navy", "white", "firebrick3"))(100)是把3种颜色切割成渐变的100种颜色，热图设置颜色都是这样设置。

##基因轨迹图
gs = head(gene_to_cluster)
plot_cell_trajectory(sc_cds,markers=gs,
                     use_color_gradient=T)

##基因拟时序点图，横坐标是按照pseudotime排好顺序的
plot_genes_in_pseudotime(sc_cds[gs,],
                         color_by = "celltype",
                         nrow= 3, #6个基因所以排了3行，数量有变化时要改
                         ncol = NULL )