（单细胞-SingleCell）拟时序分析-细胞轨迹追踪_r语言单细胞拟时序分析

在细胞轨迹分析中，仅需要准备三个数据即可

1. 原始count数据
2. 样本的meta信息
3. 基因信息

准备三个数据的主要原因是将seurat对象装换为monocle对象

创建好monocle对象后仅需要走差异基因流程和排序流程就可以进行细胞轨迹追踪

# clear
if(T){
  rm(list=ls())
  setwd("C:\\Users\\yuansh\\Desktop")
  dir()
}

# options
if(T){
  options(stringsAsFactors = F)
  options(as.is = T)
}

# packages
if(T){
  library(stringr)
  library(magrittr)
  library(ggplot2)
  library(Seurat)
  library(devtools)
  library(clustree)
  library(tidyverse)
  library(gridExtra)
  library(ggridges)
  library(ggplot2)
  library(ggExtra)
  library(monocle)
  
}

# myfunctions 
if(T){
  myhead = function(x){
    if(dim(x)[2]>5){return(x[1:5,1:5])}
    else(head(x))
  }
  myscale = function(x, scale_factor){
    (x / sum(x)) *scale_factor
  }
}


# set work-dir
setwd('C:\\Users\\yuansh\\Desktop\\DLTumorCell\\code')

# main 
# 导入数据并预处理
load("../processfile/S03_Merged_main_filtered_with_neo_osi.RData")
sce = main_tiss_filtered1
rm(main_tiss_filtered1)
if(T){
  sce@meta.data$sample_name <- as.character(sce@meta.data$sample_name)
  sample_name <- as.character(sce@meta.data$sample_name)
  # Make table 
  tab.1 <- table(sce@meta.data$sample_name) 
  # Which samples have less than 10 cells 
  samples.keep <- names(which(tab.1 > 10))
  metadata_keep <- filter(sce@meta.data, sample_name %in% samples.keep)
  # Subset Seurat object 
  sce <- subset(sce, cells=as.character(metadata_keep$cell_id))
  sce
  table(sce@meta.data$sample_name)
  table(sce@meta.data$patient_id)
  #save(sce, file=paste(dir,"S03_subset_preprocessed.RData", sep=""))
}# 剔除组织样本中少于10个细胞的组织
sce
# 导入注释信息
pre.label = read.csv("../processfile/train-predict.csv",row.names = 1)

# 这里注意一下，轨迹追中只需要上皮细胞
use.cells  <- row.names(pre.label)
sce <-subset(sce, cells=use.cells)  # 从原始数据中获取非免疫细胞表达矩阵
pre.label = pre.label$PredictLabel %>% as.data.frame() %>% set_rownames(rownames(pre.label))
sce = AddMetaData(sce,metadata = pre.label,col.name = 'prelabel')

# 轨迹追踪准备数据
# 1.原始count数据
count = sce@assays$RNA@counts %>% as.matrix()
# 2.基因信息
gene_ann <- data.frame(
  gene_short_name = row.names(count), 
  row.names = row.names(count)
)
# 3.样本注释信息
sample_ann = data.frame(
  group = sce@meta.data$prelabel,
  row.names = rownames(sce@meta.data)
)
# 初始化配置文件
pd <- new("AnnotatedDataFrame",
          data=sample_ann)
fd <- new("AnnotatedDataFrame",
          data=gene_ann)
# 穿件对象
cds <- newCellDataSet(
  count, 
  phenoData = pd,
  featureData =fd,
  expressionFamily = negbinomial.size(),
  lowerDetectionLimit=1)
cds

# 初始化对象信息
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

# 识别差异基因
# 注意这一步需要很长的时间
if(F){
  disp_table <- dispersionTable(cds)
  unsup_clustering_genes <- subset(disp_table, 
                                   mean_expression >= 0.1)
  cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id)
  dim(cds)
  diff_test_res <- differentialGeneTest(cds,fullModelFormulaStr = "~group")
  # 哪怕仅仅是65个单细胞，monocle的这个differentialGeneTest函数运行也不快。
  ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
  save(ordering_genes,file = '../processfile/ordering_genes_by_Biological_Condition_high.Rdata')
}
#导入轨迹分析需要的差异基因
# load(file = 'ordering_genes_by_Biological_Condition_high.Rdata')



cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds) 
# 然后降维
cds <- reduceDimension(cds, max_components = 2,
                       method = 'DDRTree')
# 降维是为了更好的展示数据。
# 降维有很多种方法, 不同方法的最后展示的图都不太一样, 其中“DDRTree”是Monocle2使用的默认方法
# 接着对细胞进行排序
cds <- orderCells(cds)
## 最后两个可视化函数 
plot_cell_trajectory(cds, color_by = "prelabel")  
ggsave('../plot/normal-tumor-trace.pdf',width = 8,height = 8)
# 可以很明显看到细胞的发育轨迹 
# 还有几个其它可视化函数，我们明天介绍
plot_cell_trajectory(cds, color_by = "State")
plot_cell_trajectory(cds, color_by = "Pseudotime")
plot_cell_trajectory(cds, color_by = "State") +
  facet_wrap(~State, nrow = 1)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

流程非常非常简单，感觉没什么好说的如果有不懂的可以通过以下的方式联系我

Best Regards,
Yuan.SH
---------------------------------------
School of Basic Medical Sciences,
Fujian Medical University,
Fuzhou, Fujian, China.
please contact with me via the following ways:
(a) e-mail :yuansh3354@163.com
1
2
3
4
5
6
7
8