单细胞多组学分析

单细胞转录组分析流程

单细胞技术的核心：实现单个细胞精度的分辨度

Bulk-RNA测序关注的重点是基因

单细胞技术首先关注细胞状态的差异，再关注基因表达的差异

在单细胞 RNA 测序（scRNA-seq）分析中，nFeature 和 nCount 是两个常见的质量控制（QC）指标，分别表示：

  1. nFeature（Feature Count 或 Gene Count）
    • 这是细胞中检测到的基因（或特征，通常是基因）的数量。在每个细胞中，这个值代表在这个细胞中有表达（即检测到转录本）的基因数量。如果某个细胞中基因表达的数量过少，可能表明该细胞的捕获或测序质量较低。
  2. nCount（UMI Count 或 Total Count）
    • 这是一个细胞中所有转录本的总数，通常通过唯一分子标识（UMI，Unique Molecular Identifier）来计数。UMI 可以帮助消除 PCR 扩增的重复性，以更准确地估计在每个细胞中测序到的总转录本数。

应用：

  • 在数据质量控制过程中，通常会根据这两个值过滤掉质量较差的细胞。例如，具有极低 nFeature 或 nCount 的细胞可能表示有技术噪音，反之，极高的值可能是双细胞（doublet）或其他异常情况。

Seurat 对象构建及降维聚类

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())


## install.packages('remotes')
## remotes::install_version("Seurat", version = "3.2.0")



install.packages('Seurat')
getwd()

setwd("/home/oelab/RunData2/heisenberg/Study_R/20240913/")
###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
#### 清除环境，就是把环境清空了
rm(list=ls())

## 读取10x的数据
?Read10X 
scRNA.counts=Read10X(data.dir = "/home/oelab/RunData2/heisenberg/Study_R/20240913/BC21/")
class(scRNA.counts)

### 创建Seurat对象
?CreateSeuratObject

### Seurat对象的创建主要是为了后续的分析
### min.cells 这个基因至少在3个细胞中表达
scRNA = CreateSeuratObject(scRNA.counts ,
                           min.cells = 3, project="sampl.1", 
                           min.features = 40)

view(scRNA)
View(scRNA)


### s3 data.frame  list  character matirx 
### s4 层级结构

## 第一种s4对象的提取方法：点击白框

gene=scRNA@assays[["RNA"]]@data@Dimnames[[1]]
bc=scRNA@assays[["RNA"]]@data@Dimnames[[2]]
bc[1:10]
gene[1:100]


phe=scRNA@meta.data
count=scRNA@assays[["RNA"]]@counts
z=scRNA@assays[["RNA"]]@counts@Dimnames[[2]]

##### 第二提取s4对象的方法：@ $交替使用

x1=scRNA@assays$RNA@data

count=scRNA@assays$RNA@counts
scRNA@meta.data


# 查看样本的细胞数量
table(scRNA@meta.data$orig.ident)        

## 计算质控指标

# 计算细胞中线粒体基因比例
?PercentageFeatureSet
x=scRNA[[]]
x=PercentageFeatureSet(scRNA, pattern = "^MT-")
scRNA[["percent.mt"]] <-  PercentageFeatureSet(scRNA, pattern = "^MT-")
scRNA@meta.data$per.mt=PercentageFeatureSet(scRNA, pattern = "^MT-")

# 计算红细胞比例
HB.genes <-  c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
HB_m <-  match(HB.genes, rownames(scRNA@assays$RNA)) 
HB.genes <-  rownames(scRNA@assays$RNA)[HB_m] 
HB.genes <-  HB.genes[!is.na(HB.genes)] 
scRNA[["percent.HB"]] <- PercentageFeatureSet(scRNA, features=HB.genes) 
#head(scRNA@meta.data)


a=c("a1","a2","a3")
b=c("b1","a2","m1","m2","a1")

## 匹配函数
match(a,b)


#### Feature、count、线粒体基因、红细胞基因占比可视化。
violin <-  VlnPlot(scRNA,
                     features = c("nFeature_RNA", "nCount_RNA",
                                  "percent.mt","percent.HB"), 
                     
                     pt.size = 0.01, #不需要显示点，可以设置pt.size = 0
                     ncol = 4) 
# ncol 表示一行放几个图
# pt.size 表示点的大小
violin

VlnPlot(scRNA,
        features = c("nFeature_RNA", "nCount_RNA",
                     "percent.mt","percent.HB"), 
        
        pt.size = 0.01, #不需要显示点，可以设置pt.size = 0
        ncol = 4)

?VlnPlot
###把图片画到画板上面
violin
#####以后保存图片都手动保存 不要用代码保存了。
#####你做到保存图片的时候喊师傅  我们语音聊天 我教你怎么保存


### dpi 表示图的大小 plot表示保存哪张图
ggsave("vlnplot_before_qc.pdf",dpi = 300, plot = violin, width = 12, height = 6) 
ggsave("vlnplot_before_qc.png", plot = violin, width = 12, height = 6)  

### 这几个指标之间的相关性。 把图画到画板上，然后手动保存
plot1=FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2=FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot3=FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.HB")
pearplot <- CombinePlots(plots = list(plot1, plot2, plot3), nrow=1, legend="none") 
plot1
####看画板
plot2
####看画板
plot3
####看画板
pearplot
####看画板
# 自己选择性保存图片，但是 pearplot的图片必须要看，因为这个是做质控用的
## 可以看到，nFeature_RNA的范围在0到8000之内，percent.mt代表线粒体含量
### 默认线粒体含量至少要小于20%，这是根据生物学知识得出的默认阈值。红细胞的数目要至少小于5%
### 至于nFeature_RNA和nCount_RNA的阈值怎么确定，这个要结合 pearplot的图来判断。我们质控的目标就是删除离异值。而且注意阈值尽可能取的宽松一下，防止后面分析想要的细胞得不到。
### 接下来从pearplot的图片来做质控---剔除离异值
## nFeature_RNA选择大于200 小于7500的 nFeature_RNA选择小于100000，percent.mt小于20，percent.HB小于5
scRNA1 <- subset(scRNA, subset = nFeature_RNA > 500 & 
                      percent.mt < 20 & percent.HB < 1 & nCount_RNA > 1000)
scRNA
scRNA1
### 在控制台中我们可以看到有500多细胞过滤了
#### 过滤完之后，就要对数据进行均一化，使用NormalizeData这个函数。
### 注意均一化是用NormalizeData，标准化是用ScaleData
?NormalizeData
scRNA1 <-  NormalizeData(scRNA1, normalization.method = "LogNormalize", 
                           scale.factor = 10000)


c1  g1 20  g2 40    60
c2  g1 15  g2 15    30


c1  g1 20/60  g2 40/60    1
c2  g1 15/30  g2 15/30    1



c2  g1 15  g2 15    30


### 好了，这一节数据加载、质控的内容就算是做完了。
### 在我们关闭rstudio之前 先把环境中运行好的数据保存一下
### 数据将保存在之前设定好的路径中。还有保存的scRNA1，不是scRNA，因为scRNA1才是过滤好的数据。
save(scRNA1,file='scRNA1.Rdata')


# 降维 20000 -> 2000
### 官方推荐是2000个高变基因，很多文章也有设置30000的，这个因自己的实验项目决定
scRNA1 <-  FindVariableFeatures(scRNA1, selection.method = "vst", 
                                  nfeatures = 2000) 
# Identify the 10 most highly variable genes，把top10的高变基因挑选出来，目的是为了作图
top10 <-  head(VariableFeatures(scRNA1), 10) 
# plot variable features with and without labels  画出来不带标签的高变基因图
plot1 <-  VariableFeaturePlot(scRNA1) 
### 把top10的基因加到图中
plot2 <-  LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) 
plot <-  CombinePlots(plots = list(plot1, plot2),legend="bottom") 
### 画图
plot 
#### 去画板看看 想保存就自己手动保存 


### 数据的中心化处理
### ScalaData 函数是数据比例伸缩

## 如果内存足够最好对所有基因进行中心化
scale.genes <-   rownames(scRNA1)
scRNA1 <-  ScaleData(scRNA1, features = scale.genes)

## 如果内存不够，可以只对高变基因进行标准化
scale.genes <-   VariableFeatures(scRNA1)
scRNA1 <-  ScaleData(scRNA1, features = scale.genes)


scRNA1@assays[["RNA"]]@layers[["counts"]][1:10,1:10]

scRNA1@assays[["RNA"]]@layers[["data"]][1:10,1:10]

scRNA1@assays[["RNA"]]@layers[["scale.data"]][1:10,1:10]


#细胞周期回归：上一步找到的高变基因，常常会包含一些细胞周期相关基因。
#它们会导致细胞聚类发生一定的偏移，即相同类型的细胞在聚类时会因为细胞周期的不同而分开。
?CaseMatch
cc.genes
CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA1))
#细胞周期评分
g2m_genes = cc.genes$g2m.genes
g2m_genes = CaseMatch(search = g2m_genes, match = VariableFeatures(scRNA1))
s_genes = cc.genes$s.genes
s_genes = CaseMatch(search = s_genes, match = VariableFeatures(scRNA1))
scRNA1 <-  CellCycleScoring(object=scRNA1,  g2m.features=g2m_genes,  s.features=s_genes)
#查看细胞周期基因对细胞聚类的影响
scRNAa <- RunPCA(scRNA1, features = c(s_genes, g2m_genes))
p <- DimPlot(scRNAa, reduction = "pca", group.by = "Phase")
p

VlnPlot(scRNAa, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB","G2M.Score","S.Score"), ncol = 6)
ggsave("cellcycle_pca.png", p, width = 8, height = 6)

#scRNAb <- ScaleData(scRNA1, vars.to.regress = c("S.Score", "G2M.Score"))


### PCA降维到了50个维度
scRNA1 <-  RunPCA(scRNA1 , features = VariableFeatures(object = scRNA1) )
plot1 <-  DimPlot(scRNA1, reduction = "pca", group.by="orig.ident") 
y1=scRNA1@reductions[["pca"]]@cell.embeddings

ElbowPlot(scRNA1, ndims=50, reduction="pca") 

pc.num=1:20
scRNA1 <-  FindNeighbors(scRNA1, dims = 1:20) 
scRNA1 <-  FindClusters(scRNA1, resolution = 1)  ## resolution
table(scRNA1$seurat_clusters)

scRNA1 <- BuildClusterTree(scRNA1)
PlotClusterTree(scRNA1)


### 可视化降维 t-sne umap
scRNA1 = RunTSNE(scRNA1, dims = 1:20)
embed_tsne <-  Embeddings(scRNA1, 'tsne')
write.csv(embed_tsne, 'embed_tsne.csv')
plot1 = DimPlot(scRNA1, reduction = "tsne") 
##画图
plot1
###label = TRUE把注释展示在图中
DimPlot(scRNA1, reduction = "tsne" ) 
###你会发现cluster都标了图中
?ggsave
ggsave("tSNE.pdf", plot = plot1, width = 8, height = 7)
##把图片保存一下

#UMAP---第二种可视化降维
scRNA1 <- RunUMAP(scRNA1, dims = 1:20)
embed_umap <-  Embeddings(scRNA1, 'umap')
write.csv(embed_umap,'embed_umap.csv') 
plot2 = DimPlot(scRNA1, reduction = "umap") 
DimPlot(scRNA1, reduction = "umap",label = T) 

plot2

多个样本数据整合

样本的整合 就是 去除批次效应，批次效应指：当细胞以不同的分组处理时，可能发生批次效应，批次效应可以由不同芯片上的细胞，不同测序泳道中的细胞或者不同时间收获的细胞组成，细胞经历的不同环境对转录组的测量或者转录组本身产生的影响，产生的影响存在于多个层面：实验中的细胞组之间、在同一实验室进行的实验之间或者来着不同实验室的数据集之间

Nature Genetics 处理样本整合的代码

链接：scCRC_continuum/snRNA_initial_clustering.R at main · winstonbecker/scCRC_continuum · GitHub

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
# Script to analyze scRNA data from scHTAN and scHuBMAP Projects
# For this script we will just make some qc plots and do simple clustering to divide the samples into groups
# WRB 2020

# Import libraries
library(dplyr)
library(Seurat)
library(ggplot2)
library(RColorBrewer)
library(patchwork)
library(ArchR)
library(viridis)
library(DoubletFinder)
set.seed(1)
# These are the analysis steps tht will be performed below, set execute steps to reflect which ones you want to do, note some are dependent on previous steps
# 1) Create seurat objects for individual 10x runs, run doblet finder and filter most likely doublets, merge into a seurat object containing all samples
# 2) QC and filtering
# 3) Add metadata
# 4) Normalize and scale data
# 5) Plot UMAPs
# 6) ID Cluster Markers
# 7) Initial cluster identification and seperation into compartments

execute_steps <- c(1,2,3,4,5,6,7)

# Define variables
sample_name <- "all_samples"
individual_qc_and_dublet_plot_location <- "./initial_clustering/doublet_analysis/"
analysis_parent_folder <- "./initial_clustering/"
setwd(analysis_parent_folder)
path_to_metadata <- "./hubmap_htan_metadata_atac_and_rna_final.csv"


# Define sets and locations of files for initial processing
scRNA_data_path <- "./data/"
scRNA_set_names <- c("all_samples_s12a", "all_samples_s12b", "all_samples_s5", "all_samples_s6", "all_samples_s7", "all_samples_s8", "all_samples_s3", "all_samples_s4", "all_samples_s9", "all_samples_s10", "all_samples_s11")
scRNA_sets <- list(c("EP007","A014-C-201","A002-C-010","A001-C-207","A001-C-124"), 
			c("B004-A-204", "B004-A-104", "A002-C-106", "A001-C-223", "A002-C-204", "A014-C-111"),
			c("A015-C-203","A015-C-204","A014-C-040","A002-C-201","A002-C-203"), 
			c("B001-A-301","B001-A-401","A015-C-008","A015-C-208","A014-C-114","A014-C-043","A001-C-202", "A001-C-119"),
			c("A001-C-108","A002-C-021","A002-C-212","A002-C-205","A014-C-101","A014-C-108", "A001-C-104"),
			c("A015-C-005","A015-C-006","A015-C-106","A002-C-114","A015-C-104","A015-C-202"), 
			c("A001-C-014","A001-C-023","A002-C-016","A002-C-024","A014-C-001","A014-C-054","A015-C-002","A015-C-010"),
			c("A015-C-109","B004-A-004","A002-C-121-R0", "A002-C-010-R0"), 
			c("A001-C-007", "A001-C-203", "A002-C-116", "A002-C-121", "A008-E-008", "A008-E-015", "A010-E-018", "A010-E-023", "A014-C-008", "A014-C-052"),
			c("A015-C-001", "A018-E-013", "A018-E-020", "B001-A-406", "B001-A-501", "B004-A-008", "EP034", "EP072B", "EP091", "A022-E-022"),
			c("CRC1_8810", "CRC2_15564", "CRC3_11773")) 

# Define functions
seurat_standard_normalize_and_scale <- function(colon, cluster, cluster_resolution){
	# colon is seurat object, 
	colon <- NormalizeData(colon, normalization.method = "LogNormalize", scale.factor = 10000)
	colon <- FindVariableFeatures(colon, selection.method = "vst", nfeatures = 2000)
	all.genes <- rownames(colon)
	colon <- ScaleData(colon, features = all.genes)
	colon <- RunPCA(colon, features = VariableFeatures(object = colon))
	if (cluster){
		colon <- FindNeighbors(colon, dims = 1:20)
		colon <- FindClusters(colon, resolution = cluster_resolution)
	}
	colon <- RunUMAP(colon, dims = 1:20)
	return(colon)
}

make_seurat_object_and_doublet_removal <- function(data_directory, project_name){
	# function for basic seurat based qc and doubletfinder based doublet removal
	colon.data <- Read10X(data.dir = data_directory)
	currentSample <- CreateSeuratObject(counts = colon.data, project = project_name, min.cells = 3, min.features = 40)
	currentSample[["percent.mt"]] <- PercentageFeatureSet(currentSample, pattern = "^MT-")

	# qc plot-pre filtering
	pdf(paste0("./qc_plots_", project_name, "_prefiltered.pdf"))
	print(VlnPlot(currentSample, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.05))
	dev.off()
	pdf(paste0("./qc_plots_", project_name, "_prefiltered_no_points.pdf"))
	print(VlnPlot(currentSample, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0))
	dev.off()

	# filter everything to 400 unique genes/cell
	currentSample <- subset(currentSample, subset = nFeature_RNA > 400 & nFeature_RNA < 4000)
	
	# Normalize and make UMAP
	currentSample <- seurat_standard_normalize_and_scale(currentSample, FALSE)

	# Run doublet finder
	nExp_poi <- round(0.08*length(currentSample@meta.data$orig.ident)*length(currentSample@meta.data$orig.ident)/10000)  ## Assuming 7.5% doublet formation rate - tailor for your dataset
	seu_colon <- doubletFinder_v3(currentSample, PCs = 1:20, pN = 0.25, pK = 0.09, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
	print(head(seu_colon@meta.data))
	
	# rename columns
	seu_colon$doublet.class <- seu_colon[[paste0("DF.classifications_0.25_0.09_",nExp_poi)]]
	seu_colon[[paste0("DF.classifications_0.25_0.09_",nExp_poi)]] <- NULL
	pann <- grep(pattern="^pANN", x=names(seu_colon@meta.data), value=TRUE)
	seu_colon$pANN <- seu_colon[[pann]]
	seu_colon[[pann]] <- NULL

	# plot pre and post doublet finder results
	pdf(paste0("./UMAP_pre_double_removal", project_name, ".pdf"))
	print(DimPlot(seu_colon, reduction = "umap", group.by = "doublet.class", cols = c("#D51F26", "#272E6A")))
	dev.off()
	seu_colon <- subset(seu_colon, subset = doublet.class != "Doublet")
	pdf(paste0("./UMAP_post_double_removal", project_name, ".pdf"))
	print(DimPlot(seu_colon, reduction = "umap", cols = c("#D51F26")))
	dev.off()

	# Remove extra stuff and return filtered Seurat object
	seu_colon <- DietSeurat(seu_colon, counts=TRUE, data=TRUE, scale.data=FALSE, assays="RNA")
	return(seu_colon)
}

seurat_qc_plots <- function(colon, sample_name){
	# Make some basic qc plots
	pdf(paste0("./seurat_nFeature_plots_", sample_name, ".pdf"), width = 40, height = 15)
	print(VlnPlot(colon, features = c("nFeature_RNA"), ncol = 1, pt.size = 0.2))
	dev.off()

	pdf(paste0("./seurat_nCount_plots_", sample_name, ".pdf"), width = 40, height = 15)
	print(VlnPlot(colon, features = c("nCount_RNA"), ncol = 1, pt.size = 0.2))
	dev.off()

	pdf(paste0("./seurat_pMT_plots_", sample_name, ".pdf"), width = 40, height = 15)
	print(VlnPlot(colon, features = c("percent.mt"), ncol = 1, pt.size = 0.2))
	dev.off()
}

##############################################################################################################################
#............................................................................................................................#
##############################################################################################################################
# 1) Create seurat objects for individual 10x runs, run doblet finder and filter most likely doublets, merge into a seurat object containing all samples
if (1 %in% execute_steps){
	setwd(individual_qc_and_dublet_plot_location)
	for (j in 1:length(scRNA_set_names)){
		samples <- scRNA_sets[[j]]
		print(paste0(scRNA_data_path, samples[1], "/outs/filtered_feature_bc_matrix/"))
		data_directory <- paste0(scRNA_data_path, samples[1], "/outs/filtered_feature_bc_matrix/")
		sample1 <- make_seurat_object_and_doublet_removal(data_directory, samples[1])
		seu_list <- c()
		for (i in 2:length(samples)){
			data_directory <- paste0(scRNA_data_path, samples[i], "/outs/filtered_feature_bc_matrix/")
			seu_list <- c(seu_list, make_seurat_object_and_doublet_removal(data_directory, samples[i]))
		}
		current_merge <- merge(sample1, y = seu_list, add.cell.ids = samples, project = scRNA_set_names[j])
		if (j==1){
			colon <- current_merge
		} else if (j>1){
			colon <- merge(colon, y = current_merge, project = "full_colon_project")
		}
	}
	setwd(analysis_parent_folder)
	colon[["percent.mt"]] <- PercentageFeatureSet(colon, pattern = "^MT-")
	saveRDS(colon, "initial_full_colon_proj_seurat.rds")
} elif (2 %in% execute_steps) {
	# if 1 is not a step, but 2 is then you need to load the initial object
	setwd(analysis_parent_folder)
	colon <- readRDS("initial_full_colon_proj_seurat.rds")
	# clean up if necessary:
	pann <- grep(pattern="^pANN", x=names(colon@meta.data), value=TRUE)
	for (i in pann){
	  colon[[i]] <- NULL
	}
	pann <- grep(pattern="^DF.classification", x=names(colon@meta.data), value=TRUE)
	for (i in pann){
	  colon[[i]] <- NULL
	}
}

##############################################################################################################################
#............................................................................................................................#
##############################################################################################################################
# 2) QC
if (2 %in% execute_steps){
	# create and set working directory to save qc plots
	if (!dir.exists(paste0(analysis_parent_folder, "all_samples_qc_plots"))){
		dir.create(paste0(analysis_parent_folder, "all_samples_qc_plots"))
	}
	setwd(paste0(analysis_parent_folder, "all_samples_qc_plots"))

	# make the standard seurat qc plots
	seurat_qc_plots(colon, sample_name)

	# Now subset the project (if not done already)
	colon <- subset(colon, subset = nFeature_RNA > 400 & nFeature_RNA < 4000 & percent.mt < 5 & nCount_RNA < 10000)

	# leave the qc directory
	setwd(analysis_parent_folder)
}


##############################################################################################################################
#............................................................................................................................#
##############################################################################################################################
# 3) Add metadata
if (3 %in% execute_steps){
	metadata <- read.table(path_to_metadata, header = TRUE, sep = ",", stringsAsFactors=FALSE)
	
	# remove atac column
	metadata <- metadata[,colnames(metadata)[2:28]]
	colnames(metadata) <- c("Sample", colnames(metadata)[2:27])
	metadata <- metadata[metadata$Sample != "",]

	meta_data_types <- colnames(metadata)
	for (i in 2:length(meta_data_types)){
		identities <- colon[['orig.ident']]
		for (j in 1:length(metadata$Sample)){
			identities[identities==metadata$Sample[j]] <- metadata[j,meta_data_types[i]]
		}
		colon <- AddMetaData(colon, identities$orig.ident, col.name = meta_data_types[i])
	}
}

##############################################################################################################################
#............................................................................................................................#
##############################################################################################################################
# 4) Normalize and scale data
if (4 %in% execute_steps){
	colon <- seurat_standard_normalize_and_scale(colon, TRUE, 1.0)
}

##############################################################################################################################
#............................................................................................................................#
##############################################################################################################################
# 5) Plot UMAPs
if (5 %in% execute_steps){
	# Plot by clustering, sample, and disease state
	colon <- FindClusters(colon, resolution = 0.5)
	pdf(paste0("./UMAPclustering" , ".pdf"), onefile=F)
	print(DimPlot(colon, reduction = "umap", cols = paletteDiscrete(values = unique(colon@meta.data$seurat_clusters), set = "stallion", reverse = FALSE)))
	dev.off()

	pdf(paste0("./UMAP_samples.pdf"), width = 12, onefile=F)
	print(DimPlot(colon, reduction = "umap", group.by = "orig.ident", cols = paletteDiscrete(values = unique(colon@meta.data$orig.ident), set = "stallion", reverse = FALSE)))
	dev.off()

	pdf(paste0("./UMAP_donor.pdf"), width = 12, onefile=F)
	print(DimPlot(colon, reduction = "umap", group.by = "Donor", cols = paletteDiscrete(values = unique(colon@meta.data$Donor), set = "stallion", reverse = FALSE)))
	dev.off()

	pdf(paste0("./UMAP_disease_state.pdf"), width = 6.5, onefile=F)
	print(DimPlot(colon, reduction = "umap", group.by = "DiseaseState",
		cols = c("#D51F26", "#D7CEC7", "#89288F", "#D7CEC7")) + theme_ArchR())
	dev.off()

	saveRDS(colon, "clustered_full_colon_proj_seurat.rds")
}


##############################################################################################################################
#............................................................................................................................#
##############################################################################################################################
# 6) ID Cluster Markers
if (6 %in% execute_steps){
	# Find cluster specific markers
	colon.markers <- FindAllMarkers(colon, only.pos = TRUE, min.pct = 0.2, logfc.threshold = 0.25, max.cells.per.ident = 250)
	top10 <- colon.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
}


##############################################################################################################################
#............................................................................................................................#
##############################################################################################################################
# 7) Initial cluster identification
if (7 %in% execute_steps){
	# Rough initial cluster identification
	new.cluster.ids <- c("Epithelial", #0
		"Epithelial",
		"Epithelial",
		"Epithelial",
		"Epithelial",
		"Epithelial",#5
		"Epithelial",
		"Epithelial",
		"Immune",
		"Stromal",
		"Immune",#10
		"Epithelial",
		"Stromal",
		"Epithelial",
		"Immune",
		"Epithelial",#15
		"Stromal",
		"Stromal",
		"Immune",
		"Tuft",
		"Epithelial",#20
		"Immune",
		"Epithelial",
		"Enteroendocrine",
		"Stromal",
		"Immune"#25
		)

	identities <- as.character(colon@meta.data$seurat_clusters)
	for (i in 0:length(new.cluster.ids)){
		identities[identities==as.character(i)] <- new.cluster.ids[i+1]
	}
	colon <- AddMetaData(colon, identities, col.name = "CellTypeInitial")
	pdf(paste0("./UMAP_cell_type_initial.pdf"), width = 12)
	AugmentPlot(DimPlot(colon, reduction = "umap", group.by = "CellTypeInitial"))
	dev.off()

	pdf(paste0("./UMAP_cell_type_initial_v2.pdf"), width = 12)
	plot = DimPlot(colon, reduction = "umap", group.by = "CellTypeInitial")
	plot = LabelClusters(plot = plot, id = "CellTypeInitial")
	AugmentPlot(plot)
	dev.off()

	# Save ids of immune, epithelial, and stromal cells
	immune_clusters <- c(8,10,14,18,21,25)
	stromal_clusters <- c(9,12,16,17,24)
	epithelial_clusters <- c(0:7,11,13,15,19,20,22,23)
	epithelial_specialized_clusters <- c(19,23)
	epithelial_nonspecialized_clusters <- c(0:7,12,14,20,21,25)

	# Subset to make immune, stromal, and epithelial projects
	colon_immune <- DietSeurat(subset(colon, subset = seurat_clusters %in% immune_clusters))
	colon_stromal <- DietSeurat(subset(colon, subset = seurat_clusters %in% stromal_clusters))
	colon_epitehlial <- DietSeurat(subset(colon, subset = seurat_clusters %in% epithelial_clusters))
	colon_specialized_epitehlial <- DietSeurat(subset(colon, subset = seurat_clusters %in% epithelial_specialized_clusters))
	colon_nonspecialized_epitehlial <- DietSeurat(subset(colon, subset = seurat_clusters %in% epithelial_nonspecialized_clusters))

	saveRDS(colon_immune, file = "colon_immune_all_samples_initial.rds")
	saveRDS(colon_stromal, file = "colon_stromal_all_samples_initial.rds")
	saveRDS(colon_epitehlial, file = "colon_epithelial_all_samples_initial.rds")
	saveRDS(colon_specialized_epitehlial, file = "colon_specialized_epitehlial_all_samples_initial.rds")
	saveRDS(colon_nonspecialized_epitehlial, file = "colon_nonspecialized_epitehlial_all_samples_initial.rds")
}

两种归一化算法

SCTransfrom

SCTransform 是对单细胞RNA测序数据进行归一化和方差稳定化，它通常与基于锚定点的整合方法结合使用，以实现不同数据集之间的整合，尤其是在跨样本、跨实验甚至跨物种的分析中。

基于锚定点的方法，最常见于 Seurat 的 FindIntegrationAnchors() 和 IntegrateData() 函数中，是用来在不同数据集之间找到共同的“锚定点” —— 这些点代表了不同数据集中相似的细胞或基因表达模式。

eg：

Single-cell protein activity analysis identifies recurrence-associated renal tumor macrophages

1
2
3
4
5
6
7
8
s2 <- subset(s2, subset = percent.mt < 10 & nCount_RNA > 1000 & nCount_RNA < 15000)
s2 <- SCTransform(s2,return.only.var.genes = F,verbose = T,conserve.memory = T)
## SCTransform 包含：归一化处理、找高变基因以及数据放缩

s2 <- RunPCA(s2, features = VariableFeatures(object = s2))
s2 <- RunUMAP(s2, dims = 1:30, verbose = FALSE,umap.method="umap-learn",metric="correlation")
s2 <- FindNeighbors(s2, dims = 1:30, verbose = FALSE)
s2 <- FindClusters(s2, resolution=seq(0.01,1,by=0.01), verbose = FALSE,algorithm=1)

Harmony

Harmony的优势：

  1. 整合数据的同时对稀有细胞的敏感性较好
  2. 省内存
  3. 适用于更复杂的单细胞分析实验设计，可以比较来自不同供体、组织和技术平台的细胞

eg：

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
## install.packages("devtools")
library(devtools)
## install_github("immunogenomics/harmony")
library(harmony)
#BiocManager::install("SingleCellExperiment")


scRNA.11=Read10X("BC2/")
scRNA.3=Read10X("BC21/")
scRNA.11 = CreateSeuratObject(scRNA.11 ,project="sample_11",min.cells = 3, min.features = 200)
scRNA.3 = CreateSeuratObject(scRNA.3 ,project="sample_3",min.cells = 3, min.features = 200)

scRNA_harmony <- merge(x=scRNA.11, y=c(scRNA.3 ))

scRNA_harmony <- NormalizeData(scRNA_harmony) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA(verbose=FALSE)



scRNA_harmony <- RunHarmony(scRNA_harmony, group.by.vars = "orig.ident")


#system.time({scRNA_harmony <-  RunHarmony(scRNA_harmony, group.by.vars = "orig.ident")})
###问题 harmony之后的数据在哪里？


###一定要指定harmony
scRNA_harmony <-  FindNeighbors(scRNA_harmony, reduction = "harmony", dims = 1:15) %>% FindClusters(resolution = 0.5)

scRNA_harmony <-  RunUMAP(scRNA_harmony, reduction = "harmony", dims = 1:15)
?DimPlot
plot1 =DimPlot(scRNA_harmony, reduction = "umap",label = T) 
plot2 = DimPlot(scRNA_harmony, reduction = "umap", group.by='orig.ident') 
#combinate
plotc <- plot1+plot2
plotc

现阶段单细胞流程分析总结

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
make_seurat_pp=function(x,y){
  
  scRNA.counts=Read10X(x)
  
  scRNA = CreateSeuratObject(scRNA.counts ,
                             min.cells = 3,project=y, 
                             min.features = 40)
  
  scRNA1 <- NormalizeData(scRNA, normalization.method = "LogNormalize", 
                             scale.factor = 10000)
  
  
  scRNA1 <- FindVariableFeatures(scRNA1, selection.method = "vst", 
                                    nfeatures = 2000) 
  
  
  scale.genes <- rownames(scRNA1)
  scRNA1 <- ScaleData(scRNA1, features = scale.genes)
  
  scRNA1 <-  RunPCA(scRNA1 ) 
  
  pc.num=1:20
  scRNA1 <- FindNeighbors(scRNA1, dims = pc.num) 
  scRNA1 <- FindClusters(scRNA1, resolution = 1)
  scRNA1 <- RunUMAP(scRNA1, dims = pc.num)
  
  return(scRNA1)
  
}

四种数据整合方式的对比

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
## 数据的整合

##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())

setwd("/home/oelab/RunData2/heisenberg/Study_R/20240913/")
###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)

x=list.files()
"abc"


## 内容
dir = c('./BC2/', "./BC21/")

## 起的名字
names(dir) = c('s1',  's2')      

?Read10X

## Read10X可以把两个样本读取到一个对象中
counts <- Read10X(data.dir =dir)
scRNA1 = CreateSeuratObject(counts, min.cells = 3, min.features = 200)


scRNA1[["percent.mt"]] <- PercentageFeatureSet(scRNA1, pattern = "^MT-")



VlnPlot(scRNA1,
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt" ), 
        
        pt.size = 0.01, #不需要显示点，可以设置pt.size = 0
        ncol = 3) 

scRNA2 <- subset(scRNA1, subset = nFeature_RNA > 500 & 
                      percent.mt < 20 &   nCount_RNA > 1000)



scRNA3=SplitObject(scRNA2,split.by = "orig.ident")

ifnb=scRNA2
ifnb[["RNA"]] <- split(ifnb[["RNA"]], f = ifnb$orig.ident)
ifnb


##############################################################################
##############################################################################
##############################################################################

### 如果不用整合的话，分析的流程
##Perform analysis without integration
# run standard anlaysis workflow
ifnb <-  NormalizeData(ifnb)
ifnb <-  FindVariableFeatures(ifnb)
ifnb <-  ScaleData(ifnb)
ifnb <-  RunPCA(ifnb)
ifnb <-  FindNeighbors(ifnb, dims = 1:30, reduction = "pca")
ifnb <-  FindClusters(ifnb, resolution = 1, cluster.name = "unintegrated_clusters")
ifnb <-  RunUMAP(ifnb, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")
DimPlot(ifnb, reduction = "umap.unintegrated", group.by = c("orig.ident", "seurat_clusters"))

##############################################################################
##############################################################################
##############################################################################

### 常规整合 cca
## Perform integration
ifnb <-  IntegrateLayers(object = ifnb, method = CCAIntegration, orig.reduction = "pca", new.reduction = "integrated.cca",
                           verbose =T)

# re-join layers after integration
ifnb[["RNA"]] <-  JoinLayers(ifnb[["RNA"]])

ifnb <-  FindNeighbors(ifnb, reduction = "integrated.cca", dims = 1:30)
ifnb <-  FindClusters(ifnb, resolution = 1)

ifnb <-  RunUMAP(ifnb, dims = 1:30, reduction = "integrated.cca")

# Visualization
DimPlot(ifnb, reduction = "umap", group.by = c("orig.ident","seurat_clusters"))



##############################################################################
##############################################################################
##############################################################################
##############################################################################
##############################################################################
##############################################################################

### 用SCTransform整合
## Perform integration with SCTransform-normalized datasets
ifnb <-  SCTransform(ifnb)

### SCTransform包含：NormalizeData、FindVariableFeatures、ScaleData，可以把这三个函数直接省略掉

ifnb <-  RunPCA(ifnb)
ifnb <-  RunUMAP(ifnb, dims = 1:30)
DimPlot(ifnb, reduction = "umap", group.by = c("orig.ident", "seurat_annotations"))

# integrate datasets
ifnb <-  IntegrateLayers(object = ifnb, method = CCAIntegration, normalization.method = "SCT", verbose = F)
ifnb <-  FindNeighbors(ifnb, reduction = "integrated.dr", dims = 1:30)
ifnb <-  FindClusters(ifnb, resolution = 0.6)
ifnb <-  RunUMAP(ifnb, dims = 1:30, reduction = "integrated.dr")
DimPlot(ifnb, reduction = "umap", group.by = c("orig.ident", "seurat_annotations"))

###########################################################################################


#####################################################################


### 用Harmony整合
## install.packages("devtools")
library(devtools)
## install_github("immunogenomics/harmony")
library(harmony)
#BiocManager::install("SingleCellExperiment")


scRNA.11=Read10X("BC2/")
scRNA.3=Read10X("BC21/")
scRNA.11 = CreateSeuratObject(scRNA.11 ,project="sample_11",min.cells = 3, min.features = 200)
scRNA.3 = CreateSeuratObject(scRNA.3 ,project="sample_3",min.cells = 3, min.features = 200)

scRNA_harmony <- merge(x=scRNA.11, y=c(scRNA.3 ))

scRNA_harmony <- NormalizeData(scRNA_harmony) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA(verbose=FALSE)



scRNA_harmony <- RunHarmony(scRNA_harmony, group.by.vars = "orig.ident")


#system.time({scRNA_harmony <-  RunHarmony(scRNA_harmony, group.by.vars = "orig.ident")})
###问题 harmony之后的数据在哪里？


###一定要指定harmony
scRNA_harmony <-  FindNeighbors(scRNA_harmony, reduction = "harmony", dims = 1:15) %>% FindClusters(resolution = 0.5)

scRNA_harmony <-  RunUMAP(scRNA_harmony, reduction = "harmony", dims = 1:15)
?DimPlot
plot1 =DimPlot(scRNA_harmony, reduction = "umap",label = T) 
plot2 = DimPlot(scRNA_harmony, reduction = "umap", group.by='orig.ident') 
#combinate
plotc <- plot1+plot2
plotc


################################################################################################################################

双细胞的去除

单个样本的双细胞的去除

双细胞的原理：

  1. 从现有单细胞数据中人为产生doublets
  2. 将人工产生的doublets与真实细胞混合在一起
  3. 用PCA降维或者用PCA距离矩阵寻找每个单元的artificial k最近邻居（pANN）的比例
  4. 根据预期的doublets数量进行排序和计算pANN阈值

eg：

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
##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())

setwd("/home/oelab/RunData2/heisenberg/Study_R/")
###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
#remotes::install_github('chris-mcginnis-ucsf/DoubletFinder')
library(DoubletFinder)


kidney.data=Read10X("BC10/")
## Pre-process Seurat object (standard) --------------------------------------------------------------------------------------
seu_kidney <-  CreateSeuratObject(kidney.data)
seu_kidney <-  NormalizeData(seu_kidney)
seu_kidney <-  FindVariableFeatures(seu_kidney, selection.method = "vst", nfeatures = 2000)
seu_kidney <-  ScaleData(seu_kidney)
seu_kidney <-  RunPCA(seu_kidney)
seu_kidney <-  RunUMAP(seu_kidney, dims = 1:10)


### 找pk值
## pK 控制的是每个细胞的邻域大小，即在计算双重细胞评分时，用于评估每个细胞周围的细胞数量
## pK Identification (no ground-truth) ---------------------------------------------------------------------------------------
sweep.res.list_kidney <- paramSweep(seu_kidney, PCs = 1:10, sct = FALSE)
sweep.stats_kidney <- summarizeSweep(sweep.res.list_kidney, GT = FALSE)
bcmvn_kidney <- find.pK(sweep.stats_kidney)

## Assuming 7.5% doublet formation rate - tailor for your dataset
nExp_poi <- round(0.08*length(seu_kidney@meta.data$orig.ident)*length(seu_kidney@meta.data$orig.ident)/10000)  
#DoubletRate1 = ncol(seu_kidney)*8*1e-6 

## Run DoubletFinder with varying classification stringencies ----------------------------------------------------------------
seu_kidney <- doubletFinder(seu_kidney, PCs = 1:20, pN = 0.25, pK = 0.09, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
table(seu_kidney@meta.data$DF.classifications_0.25_0.09_2445) 


seu_colon=seu_kidney
# rename columns
seu_colon$doublet.class <- seu_colon[[paste0("DF.classifications_0.25_0.09_",nExp_poi)]]
seu_colon[[paste0("DF.classifications_0.25_0.09_",nExp_poi)]] <- NULL
pann <- grep(pattern="^pANN", x=names(seu_colon@meta.data), value=TRUE)
seu_colon$pANN <- seu_colon[[pann]]
seu_colon[[pann]] <- NULL

seu_colon1 <- subset(seu_colon, subset = doublet.class == "Singlet")
seu_colon2 <- subset(seu_colon, subset = doublet.class != "Doublet")


# Remove extra stuff and return filtered Seurat object
seu_colon3 <- DietSeurat(seu_colon1, counts=TRUE, data=TRUE, scale.data=FALSE, assays="RNA")

亚群注释

目标：

  1. 确定每个类群的基因标记物
  2. 使用标记物识别每个类群的细胞类型
  3. 根据细胞类型标记物来判断是否需要重新分组，或者需要合并或者拆分聚类的类群

可以查询marker基因的网站：CellMarker (xbio.top)

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())



library(dplyr)
library(Seurat)
library(ggplot2)
library(RColorBrewer)
library(patchwork)
library(viridis)
library(DoubletFinder)
set.seed(1)

setwd("/home/oelab/RunData2/heisenberg/Study_R/")

data_directory=c("BC2/" ,"BC21/")
project_name=c("sample2", "sample21")


samples <- project_name

sample1 <- make_seurat_object_and_doublet_removal(data_directory[1], samples[1])



###  多个样本合并 
seu_list <- sample1
for (i in 2:length(samples)){
  
  
  
  sc.i = make_seurat_object_and_doublet_removal(data_directory[i], samples[i])
  seu_list=merge(seu_list,sc.i)
  
}

table(seu_list$orig.ident)


scRNA_harmony=seu_list
scRNA_harmony <- NormalizeData(scRNA_harmony ) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA(verbose=FALSE)


library(harmony)
scRNA_harmony <- RunHarmony(scRNA_harmony, group.by.vars = "orig.ident")
###问题 harmony之后的数据在哪里？


###一定要指定harmony
scRNA_harmony <- FindNeighbors(scRNA_harmony, reduction = "harmony", dims = 1:20) %&gt;% FindClusters(resolution =0.5)

scRNA_harmony <- RunUMAP(scRNA_harmony, reduction = "harmony", dims = 1:20)

DimPlot(scRNA_harmony , reduction = "umap",label = T) 
DimPlot(scRNA_harmony, reduction = "umap", split.by ='orig.ident')

DimPlot(scRNA_harmony, reduction = "umap", group.by='orig.ident')
table(scRNA_harmony$orig.ident)  

save(scRNA_harmony,file = "scRNA_harmony.rdata")


##############################################################################
##############################################################################
##############################################################################
##############################################################################
##############################################################################
##############################################################################
###  

### 另外一种整合方法  ---CCA 锚定点整合

? split
ifnb=seu_list
ifnb[["RNA"]] <- split(ifnb[["RNA"]], f = ifnb$orig.ident)
ifnb

##Perform analysis without integration
# run standard anlaysis workflow
ifnb <- NormalizeData(ifnb)
ifnb <- FindVariableFeatures(ifnb)
ifnb <- ScaleData(ifnb)
ifnb <- RunPCA(ifnb)
ifnb <- FindNeighbors(ifnb, dims = 1:30, reduction = "pca")
ifnb <- FindClusters(ifnb, resolution = 1, cluster.name = "unintegrated_clusters")
ifnb <- RunUMAP(ifnb, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")
DimPlot(ifnb, reduction = "umap.unintegrated", group.by = c("orig.ident", "seurat_clusters"))

##############################################################################

## Perform integration
ifnb <- IntegrateLayers(object = ifnb, method = CCAIntegration, orig.reduction = "pca", new.reduction = "integrated.cca",
                           verbose =T)

# re-join layers after integration
ifnb[["RNA"]] <- JoinLayers(ifnb[["RNA"]])

ifnb <- FindNeighbors(ifnb, reduction = "integrated.cca", dims = 1:30)
ifnb <- FindClusters(ifnb, resolution = 1)

ifnb <- RunUMAP(ifnb, dims = 1:30, reduction = "integrated.cca")
# Visualization
DimPlot(ifnb, reduction = "umap", group.by = c("orig.ident","seurat_clusters"))



############################################################
############################################################
############################################################


### 亚群注释


########################################################
### 方法一

### JoinLayers可以把多个样本进行合并
scRNA_harmony=JoinLayers(scRNA_harmony)

### 下载MAST包
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
BiocManager::install("MAST")
library(MAST)


### 如果找marker基因过慢，可以选择抽样
table(scRNA_harmony@meta.data$seurat_clusters)

sc.s=sample(colnames(scRNA_harmony),)

### subset 可以筛选Seurat对象的样本，downsample表示从每个分组（通常是每个细胞群或簇）中随机选择最多 100 个细胞
sc.s=subset(scRNA_harmony, downsample=100)
table(sc.s@meta.data$seurat_clusters)


### 寻找marker基因
### marker基因的特征：特异性、高表达
markers <- FindAllMarkers(sc.s, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.5,test.use = "MAST")
top10 <- markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

# Macrophage
'C1QA', 'C1QB', 'C1QC', 'APOE', 'APOC1',
#### 2  C1QA  C1QB  
scRNA_harmony=RenameIdents(scRNA_harmony,"2"="Macrophage")

scRNA_harmony=RenameIdents(scRNA_harmony, 
                           "7"="Endothelial cell", 
                           "4"="Fibroblast",
                           "1"="B cell", 
                           "5"="Myeloid",
                           "3"="NK")

table(scRNA_harmony@active.ident)
DimPlot(scRNA_harmony,label = T)


########################################################
###  方法二

table(scRNA_harmony$seurat_clusters)

Idents(scRNA_harmony)="seurat_clusters"


### 制作指定基因的气泡图
DotPlot(scRNA_harmony,features =c("LYZ","NKG7","CD3E",'CD3D','VCAN', 'CD19', 'CD79A',"PECAM1" ))


### RenameIdents：对Seurat对象进行细胞注释
scRNA_harmony=RenameIdents(scRNA_harmony,"7"="Endothelial cell","4"="Fibroblast",
                           "1"="B cell" ,"0"="Myeloid","8"="Myeloid","6"="Myeloid",
                           "2"="Myeloid", "5"="Myeloid","3"="NK_T")


### 制作细胞注释图
DimPlot(scRNA_harmony,label = T)

genes_to_check <- list(Tcells = c("CD2","CD3D","CD3E","NKG7"),
                       Bcells = c("CD79A","MS4A1"),
                       SMC = c("ACTA2","CNN1","PLN"),
                       MonoMacro = c("CD14","CD68","CD163"),
                       Dendriticcells = c("CD1C","CLEC9A","CLEC10A","LILR4"),
                       Mastcells = c("TPSAB1","TPSB2","CPA3"),
                       Fibroblasts = c("PDGFRA","FN1","DCN","LUM"),
                       Pericyte = c("MYH11","RGS5","MCAM","NOTCH3"),
                       Endothelialcells = c("PECAM1","VWF","PLVAP","CD34"),
                       Neutrophils = c('CSF3R','FCGR3B','CXCL8')) 


p_all_markers <- DotPlot(scRNA_harmony, features =genes_to_check, assay = "RNA") + 
  theme(axis.text.x = element_text(angle = 30, vjust = 0.5, hjust=1),axis.title = element_blank()) 
p_all_markers

##############################################################################

### 亚群删除

### 提取 1~10 的亚群，可以把小的亚群删除
sc1=subset(scRNA_harmony,ident=c(1:10))
DimPlot(sc1)
table(scRNA_harmony$seurat_clusters)

##############################################################################


### 方法三

# 用SingleR鉴定细胞类型

## 原理：
## 1. 计算其表达谱与各参考样本表达谱之间的Spearman相关性
## 2. 按每个标签的得分定义为相关性分布的固定分位数（默认情况下为0.8）
## 3. 对所有标签重复次操作，并将得分最高的标签作为该单细胞的注释

##############################################################################

## 文章中的代码
# use R/4.0
library(SingleR)
library(celldex)
library(Seurat)

# read in seurat object
analysis_parent_folder <- "./immune_results/"
setwd(analysis_parent_folder)
colon <- readRDS("./clustered_full_colon_immune_proj_seurat_post_doublet_filter.rds")

# set up reference, and define cell types to use
ref <- HumanPrimaryCellAtlasData()
types_to_use <- c("DC","Epithelial_cells","B_cell","Neutrophils","T_cells","Monocyte",
                  "Endothelial_cells","Neurons","Macrophage","NK_cell",
                  "BM","Platelets","Fibroblasts","Astrocyte","Myelocyte","Pre-B_cell_CD34-","Pro-B_cell_CD34+","Pro-Myelocyte")
ref <- ref[,(colData(ref)$label.main %in% types_to_use)]

# Run singleR
singler.pred <- SingleR(test = as.SingleCellExperiment(colon), ref = ref, labels = ref$label.fine)

# Add to seurat object, and then you cna plot the results if desired
colon <- AddMetaData(colon, metadata = singler.pred$labels, col.name = "SingleR.labels")

##############################################################################

BiocManager::install("SingleR")

BiocManager::install("celldex")

library(SingleR)


### HumanPrimaryCellAtlasData为人类细胞图谱计划的数据集
HumanPrimaryCellAtlasData = celldex::HumanPrimaryCellAtlasData()
save(HumanPrimaryCellAtlasData, file = "HumanPrimaryCellAtlasData.rdata")

### BlueprintEncodeData
BlueprintEncodeData = celldex::BlueprintEncodeData()
save(BlueprintEncodeData, file = "BlueprintEncodeData.rdata")

DatabaseImmuneCellExpressionData = celldex::DatabaseImmuneCellExpressionData()
save(DatabaseImmuneCellExpressionData, file = "DatabaseImmuneCellExpressionData.rdata")

ImmGenData = celldex::ImmGenData()
save(ImmGenData,file = "ImmGenData.rdata")

MonacoImmuneData = celldex::MonacoImmuneData()
save(MonacoImmuneData,file = "MonacoImmuneData.rdata")


MouseRNAseqData = celldex::MouseRNAseqData()
save(MouseRNAseqData,file = "MouseRNAseqData.rdata")

NovershternHematopoieticData = celldex::NovershternHematopoieticData()
save(NovershternHematopoieticData,file = "NovershternHematopoieticData.rdata")


### 下载好数据库后，把ref_Human_all.Rdata加载到环境中，这样算是对数据库的加载，就可以按照singler的算法来对细胞亚群进行定义了。
load("F:/ref.data/ref_Human_all.RData")

### 可以看到在环境中多了一个叫ref_Human_all的文件 大小为113mb  这个就是数据库
#### 然后把环境中的ref_Human_all赋值与refdata
refdata <- HumanPrimaryCellAtlasData

### 把RNA的转录表达数据提取
?LayerData
testdata <-  LayerData(scRNA_harmony, assay = "RNA", layer = "data")

### 把scRNA数据中的seurat_clusters提取出来，注意这里是因子类型的
clusters <-  scRNA_harmony@meta.data$seurat_clusters

table(refdata$label.main)

### 开始用singler分析
cellpred <- SingleR(test = testdata, ref = refdata, labels = refdata$label.main, 
                       method = "cluster", clusters = clusters, 
                       assay.type.test = "logcounts", assay.type.ref = "logcounts")


### 制作细胞类型的注释文件
celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels)

### 保存一下
write.csv(celltype,"celltype_singleR.csv",row.names = FALSE)

##把singler的注释写到metadata中 有两种方法

### 方法一
scRNA_harmony@meta.data$celltype ="NA"
for(i in 1:nrow(celltype)){
  scRNA_harmony@meta.data[which(scRNA_harmony@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]
}


### 第一次循环  
i = 1
="Macrophage"
celltype$celltype[1]

[which(scRNA_harmony@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype']

celltype$ClusterID[1]
which(scRNA_harmony@meta.data$seurat_clusters == celltype$ClusterID[i])

which(scRNA_harmony@meta.data$seurat_clusters ==  "0")
a=c("a1","a2","a3")
b=c("b1","a2","m1","m2","a1")
which(b=="a1")


### 因为把singler的注释加载到metadata中时候，命名的名字叫celltype，所以画图时候，group.by="celltype"
DimPlot(scRNA_harmony, group.by="celltype", label=T, label.size=5)

### 方法二：
celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels ) 
scRNA_harmony@meta.data$singleR=celltype[match(clusters,celltype$ClusterID),'celltype']

### 因为把singler的注释加载到metadata中时候，命名的名字叫singleR，所以画图时候，group.by="singleR"
DimPlot(scRNA_harmony, group.by="singleR", label=T, label.size=5)
### 可以看到两种方法得到的结果都是一样的，但是我比较喜欢第二种方法

## 鉴定结果展示
DimPlot(scRNA_harmony, group.by="celltype", label=T, label.size=5, reduction='umap')

### 方法三：Azimuth
### 细胞预测
### https://satijalab.org/azimuth/

亚群细分

就是在亚群注释的基础上，单独分析亚群注释上几个cluster

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
## 亚群细分

##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())


###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)

load("/home/oelab/RunData2/heisenberg/scRNA_harmony.rdata")
setwd("/home/oelab/RunData2/heisenberg/AllData/")


### 多个样本的整合要用到JoinLayers这个函数
scRNA_harmony=JoinLayers(scRNA_harmony)

# PTPRC表示免疫细胞

scRNA1 = scRNA_harmony
Idents(scRNA1) = "seurat_clusters"

## 针对PTPRC的，提取出来的cluster
DotPlot(scRNA1, features = "PTPRC", group.by ="seurat_clusters" )

## 针对EPCAM的，提取出来的cluster
DotPlot(scRNA1, features = "EPCAM", group.by ="seurat_clusters" )

## 将PTPRC，提取出来的cluster用一个对象来接收
sc.immu=subset(scRNA1, idents = c(0,2,3,5,6,8,10))

## 绘制PTPRC的基因分布
FeaturePlot(sc.immu, features ="PTPRC")

DimPlot(scRNA1, label = T)

table(scRNA1@active.ident)
################################################################################
#sc.immu=subset(scRNA1,ident=c("immu"))
#sc.immu=subset(scRNA1,ident=c("0","2"))

## 进行harmony整合
scRNA_harmony1 <- sc.immu
scRNA_harmony1 <- NormalizeData(scRNA_harmony1) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA(verbose=FALSE)
library(harmony)
system.time({scRNA_harmony1 <- RunHarmony(scRNA_harmony1, group.by.vars = "orig.ident")})

### 一定要指定harmony
scRNA_harmony1 <- FindNeighbors(scRNA_harmony1, reduction = "harmony", dims = 1:15) %>% FindClusters(resolution = 0.5)

scRNA_harmony1 <- RunUMAP(scRNA_harmony1, reduction = "harmony", dims = 1:15)

plot1 = DimPlot(scRNA_harmony1, reduction = "umap",label = T) 
plot2 = DimPlot(scRNA_harmony1, reduction = "umap", group.by='orig.ident') 
# combinate
plotc <- plot1+plot2

## 这个结果就是针对免疫细胞的再次降维聚类 
plotc


markers <- FindAllMarkers(scRNA_harmony1, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.5,test.use = "MAST")
top10 <- markers %&>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

Marker基因展示

样本分类提取的方法

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())

###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(harmony)
library(cowplot)
library(ggrepel)
library(reshape2)

setwd("/home/oelab/RunData2/heisenberg/AllData/")
load("scRNA1.rdata")


DimPlot(scRNA1, label = T, reduction = "umap")

colnames(scRNA1@meta.data)

table(scRNA1@meta.data$orig.ident)

DimPlot(scRNA1,label = T, reduction = "umap", group.by = "orig.ident")

## ncol 就是 2*2 的展示方式
DimPlot(scRNA1,label = T, reduction = "umap", split.by  = "orig.ident", ncol = 2)


DimPlot(scRNA1,label = T, reduction = "umap")
Idents(scRNA1)="orig.ident"
DimPlot(scRNA1,label = T, reduction = "umap")

DimPlot(scRNA1,label = T, reduction = "umap", split.by  = "orig.ident" )


colnames(scRNA1@meta.data)
table(scRNA1@meta.data$orig.ident)
Idents(scRNA1)="orig.ident"
DimPlot(scRNA1,label = T)

Idents(scRNA1)="celltype"
DimPlot(scRNA1,label = T)
DimPlot(scRNA1,split.by = "orig.ident")

table(scRNA1$orig.ident)
DimPlot(scRNA1,label = T,group.by = "orig.ident")




DimPlot(scRNA1,split.by = "orig.ident")
DimPlot(scRNA1,group.by = "orig.ident")

### 参考78-79 https://github.com/UMCCR-RADIO-Lab/snRNA-seq-atlas-of-pheochromocytoma-and-paraganglioma/blob/main/snRNA_add_and_format_metadata.R

## 展示n个分组的结果

###################################
### 方法一
metadata=scRNA1@meta.data
table(scRNA1@meta.data$orig.ident)
meta.caf1 = filter(metadata, orig.ident=="CAF1")
meta.caf1$group = "CAF"
meta.caf2 = filter(metadata, orig.ident=="CAF2")
meta.caf2$group ="CAF"
meta.dap1 = filter(metadata, orig.ident=="DapiNeg1")
meta.dap1$group = "DAP"
meta.dap2 = filter(metadata, orig.ident=="DapiNeg2")
meta.dap2$group = "DAP"
meta.all = rbind(meta.caf1, meta.caf2, meta.dap1, meta.dap2)
### 参考78-79 https://github.com/UMCCR-RADIO-Lab/snRNA-seq-atlas-of-pheochromocytoma-and-paraganglioma/blob/main/snRNA_add_and_format_metadata.R

scRNA1=AddMetaData(scRNA1, metadata = meta.all)

DimPlot(scRNA1, split.by = "group", group.by = "celltype")
VlnPlot(scRNA1, features = "Rbp1", split.by = "group")

##################################################################################

### 方法二
Idents(scRNA1)="orig.ident"
table(scRNA1$orig.ident)

scRNA1=RenameIdents(scRNA1, "CAF1"="CAF", "CAF2"="CAF",
                    "DapiNeg1"="DapiNeg", "DapiNeg2"="DapiNeg")


table(scRNA1@active.ident)

scRNA1@meta.data$group.2=scRNA1@active.ident
DimPlot(scRNA1,split.by = "group.2",group.by = "celltype")

#####################################################################################
### 方法三（用的最多）

scRNA1@meta.data[which(scRNA1@meta.data$orig.ident=="CAF1"),"group.3"]="CAF"
scRNA1@meta.data[which(scRNA1@meta.data$orig.ident=="CAF2"),"group.3"]="CAF"

scRNA1@meta.data[which(scRNA1@meta.data$orig.ident=="DapiNeg1"),"group.3"]="DapiNeg"
scRNA1@meta.data[which(scRNA1@meta.data$orig.ident=="DapiNeg2"),"group.3"]="DapiNeg"

DimPlot(scRNA1,split.by = "group.3",group.by = "celltype")

####################################################################################
### 方法四

## 字符串的提取
substr("abcdef", 2, 4)

scRNA1$group.4 = substr(scRNA1@meta.data$orig.ident,1,3)

DimPlot(scRNA1, split.by = "group.4", group.by = "celltype")


####################################################################################
####################################################################################
####################################################################################


## 对于cluster的提取

DimPlot(scRNA1)
colnames(scRNA1@meta.data)
Idents(scRNA1)="seurat_clusters"
sc.1=subset(scRNA1,ident=c(0:8,10,11))
DimPlot(sc.1,label = T)


## 对于基因的提取
sc.3=subset(sc.1, Rbp1 > 0 & Gapd > 0)
table(scRNA1$orig.ident)
Idents(scRNA1) = "celltype"
DimPlot(scRNA1,label = T)

table(scRNA1@active.ident)

## 对于downsample的提取
sc.4=subset(scRNA1, downsample=50)
table(sc.4@active.ident)


scRNA2 <- subset(scRNA1, subset = nFeature_RNA > 500 & 
                      percent.mt < 20 & nCount_RNA > 1000)


##################################################
### 提取免疫细胞 是否合理？

## 不能根据参数进行细胞的提取，有可能存在双细胞
sc.cd45=subset(scRNA1, Ptprc > 0  )

## 只能根据cluster进行细胞的提取
DotPlot(scRNA1,features ="Ptprc" )
sc.cd45=subset(scRNA1,idents = c(0,5,6,9,10,13:17,21) )



Idents(scRNA1) = "celltype"
DimPlot(scRNA1, label = T)

table(scRNA1@active.ident)
Idents(scRNA1) = "celltype"
table(scRNA1@active.ident)
sc.fibr = subset(scRNA1,idents = c("Fibroblasts"))

## 对于Seurat数据集的切割
sc.fibr.t = subset(scRNA1,idents = c("Fibroblasts", "T cells"))

Marker基因绘图展示

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568

### Marker基因的展示


##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())


###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)

setwd("/home/oelab/RunData2/heisenberg/AllData/")


## 加载harmony数据
load("scRNA_harmony.rdata")

###########################################################################################################

### 绘图

pbmc3k.final=scRNA_harmony

library(ggplot2)
library(patchwork)

features <- c("LYZ", "CCL5", "IL32", "GAPDH")

## 峰峦图
RidgePlot(pbmc3k.final, features = features, ncol = 2)

## 小提琴图
## pt.size=0 就是把点给去掉
VlnPlot(pbmc3k.final, features = features, ncol = 2, pt.size=0)
?VlnPlot
## 将有表达的基因（eg：细胞组LYZ表达量大于1的）挑选出来
VlnPlot(subset(pbmc3k.final,LYZ > 0), features ="LYZ",pt.size=0) 


## 维度图
## reduction = "tsne" 这里可以tsne，也可以umap
FeaturePlot(pbmc3k.final, features = features, reduction = "tsne")

pbmc3k.final=RunTSNE(pbmc3k.final, reduction = "harmony", dims = 1:20)
FeaturePlot(pbmc3k.final, features = features, reduction = "tsne")


## 气泡图
## RotatedAxis()可以将基因旋转45°，防止重叠
DotPlot(pbmc3k.final, features = features, scale=F)
DotPlot(pbmc3k.final, features = features, scale=F) + RotatedAxis()
?DotPlot

## 热图
DoHeatmap(pbmc3k.final, features = features, size = 3)
DoHeatmap(pbmc3k.final, features = features, size = 3, slot="data")

## downsample = 100 会保持每个cluster的数量一致
DoHeatmap(subset(pbmc3k.final, downsample = 100), features = features, size = 3)
?DoHeatmap

## 查看某个基因的维度图
FeaturePlot(pbmc3k.final, features = "MS4A1")
## min.cutoff 表达值的最低值
## max.cutoff 表达值得最高值
FeaturePlot(pbmc3k.final, features = "MS4A1", min.cutoff = 1, max.cutoff = 2)
FeaturePlot(pbmc3k.final, features = c("MS4A1", "GAPDH"), min.cutoff = "q10", max.cutoff = "q90")
FeaturePlot(pbmc3k.final, features = c("MS4A1", "CD79A"), blend = TRUE)
FeaturePlot(pbmc3k.final, features = c("MS4A1", "CD79A"), split.by = "orig.ident")

### 降维图
plot1 <- DimPlot(pbmc3k.final)

## 查看两个基因是否线性相关，用FeatureScatter进行展示
plot2 <- FeatureScatter(pbmc3k.final, feature1 = "LYZ", feature2 = "CCL5")

# Combine two plots
plot1 + plot2

baseplot <- DimPlot(pbmc3k.final, reduction = "umap")

# 添加标签和标题
baseplot + labs(title = "huage PBMCs")+xlab("x1")
BiocManager::install("ggmin")
baseplot + ggmin::theme_powerpoint()


#########################################################################################################


## metadata分组


##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())

setwd("/home/oelab/RunData2/heisenberg/AllData/")

###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(harmony)
library(cowplot)

dir=c("CAFs.1/","CAFs.2/","dapi1/","dapi2/")
dir=c("CAFs.1/","CAFs.2/","dapi1/","dapi2/")
names(dir) = c('CAF1',  'CAF2', 'DapiNeg1',"DapiNeg2")      



counts <- Read10X(data.dir =dir)
scRNA1 = CreateSeuratObject(counts,min.cells = 3, min.features = 200)

scRNA1[["percent.mt"]] <- PercentageFeatureSet(scRNA1, pattern = "^mt-")



VlnPlot(scRNA1,
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt" ), 
        
        pt.size = 0.01, #不需要显示点，可以设置pt.size = 0
        ncol = 3) 

scRNA2 <- subset(scRNA1, subset = nFeature_RNA > 500 & 
                      percent.mt < 20 & nCount_RNA > 1000)



scRNA3 = SplitObject(scRNA2 , split.by = "orig.ident")

ifnb=scRNA2
ifnb[["RNA"]] <- split(ifnb[["RNA"]], f = ifnb$orig.ident)
ifnb

##Perform analysis without integration
# run standard anlaysis workflow
ifnb <- NormalizeData(ifnb)
ifnb <-  FindVariableFeatures(ifnb)
ifnb <-  ScaleData(ifnb)
ifnb <-  RunPCA(ifnb)
ifnb <-  FindNeighbors(ifnb, dims = 1:30, reduction = "pca")
ifnb <-  FindClusters(ifnb, resolution = 1, cluster.name = "unintegrated_clusters")
ifnb <-  RunUMAP(ifnb, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated")
DimPlot(ifnb, reduction = "umap.unintegrated", group.by = c("orig.ident", "seurat_clusters"))

##############################################################################
##############################################################################
##############################################################################

## Perform integration
ifnb <- IntegrateLayers(object = ifnb, method = CCAIntegration, orig.reduction = "pca", new.reduction = "integrated.cca",
                           verbose =T)

# re-join layers after integration
ifnb[["RNA"]] <- JoinLayers(ifnb[["RNA"]])

ifnb <- FindNeighbors(ifnb, reduction = "integrated.cca", dims = 1:30)
ifnb <- FindClusters(ifnb, resolution = 1)

ifnb <- RunUMAP(ifnb, dims = 1:30, reduction = "integrated.cca")
# Visualization
DimPlot(ifnb, reduction = "umap", group.by = c("orig.ident", "seurat_clusters"))

###################################################################
scRNA1=ifnb

DefaultAssay(scRNA1) <- "RNA"
##BiocManager::install("SingleR")
library(SingleR)
load("F:/ref.data/ref_Mouse_all.RData")
refdata <- ref_Mouse
###把rna的转录表达数据提取
testdata=LayerData(scRNA1, assay="RNA", layer='data')

###把scRNA数据中的seurat_clusters提取出来，注意这里是因子类型的
clusters <- scRNA1@meta.data$seurat_clusters
###开始用singler分析
cellpred <- SingleR(test = testdata, ref = refdata, labels = refdata$label.main, 
                       method = "cluster", clusters = clusters, 
                       assay.type.test = "logcounts", assay.type.ref = "logcounts")

celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels, stringsAsFactors = FALSE)
scRNA1@meta.data$celltype = "NA"
for(i in 1:nrow(celltype)){
  scRNA1@meta.data[which(scRNA1@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]
}



DimPlot(scRNA1, reduction = "umap", group.by = "celltype",label = T)
DimPlot(scRNA1, reduction = "umap", label = T)
 
getwd()
setwd("/home/oelab/RunData2/heisenberg/AllData/")
save(scRNA1,file = "scRNA1.rdata")

load("scRNA1.rdata")
#################################################################################################
##scRNA1@active.ident=scRNA1@meta.data$celltype


### 做图


my36colors <- c('#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3', '#476D87',
                  '#E95C59', '#E59CC4', '#AB3282', '#23452F', '#BD956A', '#8C549C', '#585658',
                  '#9FA3A8', '#E0D4CA', '#5F3D69', '#C5DEBA', '#58A4C3', '#E4C755', '#F7F398',
                  '#AA9A59', '#E63863', '#E39A35', '#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B',
                  '#712820', '#DCC1DD', '#CCE0F5',  '#CCC9E6', '#625D9E', '#68A180', '#3A6963',
                  '#968175'
)
DimPlot(scRNA1, group.by = "seurat_clusters", label=T, cols = my36colors, reduction = "umap")
DimPlot(scRNA1,  label = T, reduction = "umap")


DimPlot(scRNA1, group.by = "celltype", reduction = "tsne")
DimPlot(scRNA1, group.by = "celltype", reduction = "umap")
DimPlot(scRNA1)
DimPlot(scRNA1, group.by = "celltype")


colnames(scRNA1@meta.data)
Idents(scRNA1)= "celltype"
DimPlot(scRNA1)
Idents(scRNA1)="orig.ident"
DimPlot(scRNA1)

Idents(scRNA1)= "celltype"
degs <- FindAllMarkers(scRNA1, logfc.threshold = 0.5,
                          test.use = "roc", 
                          return.thresh = 0.25, 
                          min.pct = 0.3, only.pos = T) 

degs_sig <-  degs %>% 
  filter(pct.1 &gt; 0.3 &
           power &gt; 0.25) %>%
  filter(cluster != "other") %>%
  arrange(cluster, -power)  

# select degs for heatmap
degs_top50 <- degs_sig %>% 
  #  filter(cluster!="other") %>%
  group_by(cluster) %>% 
  top_n(50, power) %>%
  top_n(50, avg_diff) %>%
  arrange(cluster, -power)


avgData <- scRNA1@assays$RNA@data[degs_top50$gene,] %>% 
  apply(1, function(x){
    tapply(x, scRNA1$celltype, mean) # ExpMean
  }) %>% t


phData <- MinMax(scale(avgData), -2, 2) # z-score
rownames(phData) <- 1:nrow(phData)

##install.packages("pheatmap")
library(pheatmap)

phres <- pheatmap(
  phData, 
  color = colorRampPalette(c("darkblue", "white", "red3"))(99), #配色
  scale = "row",
  cluster_rows = F, #不按行聚类
  cluster_cols = T, #按列聚类
  clustering_method = "complete",
  show_rownames = F, #显示cluster名
  annotation_row = data.frame(cluster = degs_top50$cluster), 
)  

phData1= phData[, c("Fibroblasts","Adipocytes","Macrophages","Granulocytes","Monocytes","T cells")] 


phres <- pheatmap(
  phData1, 
  color = colorRampPalette(c("darkblue", "white", "red3"))(99), #配色
  scale = "row",
  cluster_rows = F, #不按行聚类
  cluster_cols = F, #按列聚类
  clustering_method = "complete",
  show_rownames = F, #显示cluster名
  annotation_row = data.frame(cluster = degs_top50$cluster), 
) 

phres
########################################################################################################
###################################################################################################


load("scRNA1.rdata")

Idents(scRNA1)="celltype"

scRNA1=scRNA2

VlnPlot(scRNA1,features = "Gapdh")
scRNA1@active.ident=factor(scRNA1@active.ident,levels = c( "T cells","Fibroblasts" ,
                                                           "Granulocytes", "Macrophages",   "Monocytes","Adipocytes"))
scRNA1$celltype1=scRNA1@active.ident

DimPlot(scRNA1, reduction = "umap", group.by = "celltype1",label = T)
DimPlot(scRNA1  )
VlnPlot(scRNA1,features = "Gapdh")


DotPlot(scRNA1,features = c("Cfd","S100a8","Gapdh"))



scRNA1@active.ident=factor(scRNA1@active.ident,levels = c("Adipocytes","Fibroblasts" , "Granulocytes", "Macrophages",   "Monocytes","T cells"))
table(scRNA1@active.ident)


marker=FindAllMarkers(scRNA1,logfc.threshold = 0.5)

write.csv(marker,file = "marker.csv")

top10 <- marker  %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

scRNA1 <- ScaleData(scRNA1,features = rownames(scRNA1))
DoHeatmap(scRNA1,features = top10$gene,label = F,slot = "scale.data")
table(scRNA1@active.ident)
DoHeatmap( subset(scRNA1,downsample=300) ,features = top10$gene,label = F,slot = "scale.data")


table(scRNA1$celltype)
scRNA2=scRNA1

table(scRNA1@active.ident)
scRNA1=subset(scRNA1, downsample = 150)
library(pheatmap)

colanno=scRNA1@meta.data 
colanno$barcode=rownames(colanno)
colanno=colanno%&gt;%arrange(celltype)
rownames(colanno)=colanno$barcode
colanno$barcode=NULL



colanno1=colanno[,6]
colanno1=as.data.frame(colanno1)
rownames(colanno1)=rownames(colanno)
colnames(colanno1)="celltype"

rowanno=top10
colnames(top10)
rowanno=rowanno%>%arrange(cluster)

mat4=scRNA1[["RNA"]]@scale.data[rowanno$gene,rownames(colanno1)]
mat4[mat4>=2.5]=2.5
mat4[mat4 < (-1.5)]= -1  

pheatmap(mat4,cluster_rows = F,cluster_cols = F,
         show_colnames = F,
         annotation_col = colanno1,
         gaps_row=as.numeric(cumsum(table(rowanno$cluster))[-6]),
         gaps_col=as.numeric(cumsum(table(colanno$celltype))[-6]) 
         
)



pheatmap(mat4,cluster_rows = F,cluster_cols = F,
         show_colnames = F,
         annotation_col = colanno1,
         gaps_row=as.numeric(cumsum(table(rowanno$cluster))[-6]),
         gaps_col=as.numeric(cumsum(table(colanno$celltype))[-6]),
         color = colorRampPalette(c("darkblue", "white", "red3"))(99)
         
)



pheatmap(mat4,cluster_rows = F,cluster_cols = F,
         show_colnames = F,
         annotation_col = colanno1,
         gaps_row=as.numeric(cumsum(table(rowanno$cluster))[-6]),
         gaps_col=as.numeric(cumsum(table(colanno$celltype))[-6]) ,
         filename="marker.heatmap.pdf",width=11,height = 7
)


################################################################################################


######################################################################################################
VlnPlot(scRNA1, features = top10$gene[1:3], pt.size = 0, ncol = 1)+
  scale_x_discrete("")+
  theme(
    axis.text.x.bottom = element_blank()
  )


top10 <-  marker  %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)


library(reshape2)
vln.df=as.data.frame(scRNA1[["RNA"]]@data[top10$gene,])
vln.df$gene=rownames(vln.df)
vln.df=melt(vln.df,id="gene")
colnames(vln.df)[c(2,3)]=c("barcode","exp")

anno=scRNA1@meta.data
anno$barcode=rownames(anno)


vln.df=inner_join(vln.df,anno,by="barcode")
##vln.df$gene=factor(vln.df$gene,levels = top10$gene[1:7]) #为了控制画图的基因顺序

vln.df%&>%ggplot(aes(celltype,exp))+geom_violin(aes(fill=gene),scale = "width")+
  facet_grid(vln.df$gene~.,scales = "free_y")+
  scale_fill_brewer(palette = "Set3",direction = 1)+
  scale_x_discrete("")+scale_y_continuous("")+
  theme_classic()+
  theme(
    axis.text.x.bottom = element_text(angle = 45,hjust = 1,vjust = 1),
    panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
    legend.position = "none"
  )



###########################################################################################################
## 
top10 <-  marker  %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
DotPlot(scRNA1, features = top10$gene,cols = c("lightgrey", "#FF00FF"))+
  scale_x_discrete("huage")+scale_y_discrete("")+theme_classic()+RotatedAxis()
# 
?DotPlot
bubble.df=as.matrix(scRNA1[["RNA"]]@data[ top10$gene,])
bubble.df=t(bubble.df)
bubble.df=as.data.frame(scale(bubble.df))
bubble.df$CB=rownames(bubble.df)
meta=scRNA1@meta.data[,c(1,6)]
bubble.df=merge(bubble.df,meta,by.x = "CB",by.y = 0)
bubble.df$CB=NULL

celltype_v=c()
gene_v=c()
mean_v=c()
ratio_v=c()
for (i in unique(bubble.df$celltype)) {
  bubble.df_small=bubble.df%>%filter(celltype==i)
  for (j in top10$gene) {
    exp_mean=mean(bubble.df_small[,j])
    exp_ratio=sum(bubble.df_small[,j] &gt; min(bubble.df_small[,j])) / length(bubble.df_small[,j])
    celltype_v=append(celltype_v,i)  ##Add elements to a vector.
    gene_v=append(gene_v,j)
    mean_v=append(mean_v,exp_mean)
    ratio_v=append(ratio_v,exp_ratio)
  }
}

plotdf=data.frame(
  celltype=celltype_v,
  gene=gene_v,
  exp=mean_v,
  ratio=ratio_v
)
library(RColorBrewer)
plotdf$celltype=factor(plotdf$celltype,levels = sort(unique(plotdf$celltype)))
plotdf$gene=factor(plotdf$gene,levels = rev(as.character(top10$gene[1:10])))
plotdf$exp=ifelse(plotdf$exp&gt;3,3,plotdf$exp)
plotdf%&gt;%ggplot(aes(x=celltype,y=gene,size=ratio,color=exp))+geom_point()+
  scale_x_discrete(" ")+scale_y_discrete(" ")+
  scale_color_gradientn(colours = rev(c("#FFD92F","#FEE391",brewer.pal(11, "Spectral")[7:11])))+
  scale_size_continuous(limits = c(0,1))+theme_classic()+
  theme(
    axis.text.x.bottom = element_text(hjust =1, vjust = 1, angle = 90)
  )






################################################################################
scRNA1=scRNA2

FeaturePlot(scRNA1,features = "Mzb1",reduction = "tsne",pt.size = 1,split.by = "orig.ident")
FeaturePlot(scRNA1,features = "S100a9",reduction = "tsne",pt.size = 1 )

FeaturePlot(scRNA1,features = "Rbp1",reduction = "tsne",pt.size = 1)+
  scale_x_continuous("")+scale_y_continuous("")+
  theme_bw()+ 
  theme(  
    panel.grid.major = element_blank(),panel.grid.minor = element_blank(),  
    axis.ticks = element_blank(),axis.text = element_blank(),  
    legend.position = "none",  
    plot.title = element_text(hjust = 0.5,size=14)  
  )

#### 

mat1=as.data.frame(scRNA1[["RNA"]]@data["Rbp1",])
colnames(mat1)="exp"
mat2=Embeddings(scRNA1,"tsne")
mat3=merge(mat2,mat1,by="row.names")


mat3%>%ggplot(aes(tSNE_1,tSNE_2))+geom_point(aes(color=exp))+
  scale_color_gradient(low = "grey",high ="#1B9E77")+theme_bw()+ 
  theme(  
    panel.grid.major = element_blank(),panel.grid.minor = element_blank(),  
    axis.ticks = element_blank(),axis.text = element_blank(),  
    legend.position = "none",  
    plot.title = element_text(hjust = 0.5,size=14)  
  )+
  scale_x_continuous("")+scale_y_continuous("")+labs(title= "Rbp1")




library(ggpubr)

VlnPlot( scRNA1 , "Rbp1",group.by = "orig.ident")+ 
  theme_bw()+stat_compare_means(method = 't.test')+geom_boxplot()

VlnPlot(subset(scRNA1,Rbp1 > 0 ), "Rbp1",group.by = "orig.ident")+ 
  theme_bw()+stat_compare_means(method = 't.test')+geom_boxplot()

sc.FURIN=subset(scRNA1,Rbp1 > 0 )
FURIN.exp=sc.FURIN@assays$RNA@data["Rbp1",]
FURIN.exp.m= sc.FURIN@meta.data
table(sc.FURIN$orig.ident)
FURIN.exp.m$exp=FURIN.exp
FURIN.exp.m$Cohort=factor(FURIN.exp.m$orig.ident,
                          levels = c("DapiNeg1","DapiNeg2","CAF1",'CAF2'))
write.csv(FURIN.exp.m,file = "FURIN.exp.m.csv")
colnames(FURIN.exp.m)
library(ggplot2)
library("ggsignif")
my_comparisons=list(c("DapiNeg1","CAF1"),c("DapiNeg2",'CAF2'),
                    c("DapiNeg2" ,"CAF1"))

ggplot(FURIN.exp.m,aes(Cohort,exp,fill=Cohort)) +
  geom_violin()+theme_classic()+
  theme(plot.title=element_text(size = 15,face="bold"),
        axis.text.x=element_text(size=15,face="bold",angle=25,hjust = 1,vjust = 1),
        axis.text.y=element_text(size=18,face="bold"),
        axis.title.x=element_text(size = 0,vjust = 1,face="bold"),
        axis.title.y=element_text(size = 15,face="bold"))+
  labs(y= 'Expression',x="Group")+
  geom_signif(comparisons = my_comparisons,step_increase = 0.1,
              map_signif_level = F,test = t.test,size=1,textsize = 6)+
  geom_jitter(shape=16, position=position_jitter(0.2))+
  ggtitle("Rbp1")+
  theme(plot.title = element_text(size =18,hjust = 0.5))+
  guides(fill=guide_legend(title = " "))+ ##更改legend名字
  ##如果是color分组 需要 guides(color=guide_legend(title = "Group")) 
  theme(axis.line.x=element_line(linetype=1,color="black",size=1))+
  theme(axis.line.y=element_line(linetype=1,color="black",size=1))+
  theme(legend.text = element_text(size = 15) ) ##legend字体变大

个性化绘图

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632

##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())

###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(harmony)
library(cowplot)

setwd("D:/0_SubjectData/")
load("scRNA_harmony.rdata")
DimPlot(scRNA_harmony)


############################################################
############################################################
############################################################
scRNA_harmony=JoinLayers(scRNA_harmony)
library(SingleR)

load("HumanPrimaryCellAtlasData.rdata")
refdata <- HumanPrimaryCellAtlasData
testdata = LayerData(scRNA_harmony, assay="RNA", layer='data')

### 把scRNA数据中的seurat_clusters提取出来，注意这里是因子类型的
clusters <- scRNA_harmony@meta.data$seurat_clusters

table(refdata$label.main)
### 开始用singler分析
cellpred <- SingleR(test = testdata, ref = refdata, labels = refdata$label.main, 
                       method = "cluster", clusters = clusters, 
                       assay.type.test = "logcounts", assay.type.ref = "logcounts")
### 制作细胞类型的注释文件
celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels)

scRNA_harmony@meta.data$celltype ="NA"
for(i in 1:nrow(celltype)){
  scRNA_harmony@meta.data[which(scRNA_harmony@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] &lt;- celltype$celltype[i]
}
DimPlot(scRNA_harmony)

my36colors <- c('#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3', '#476D87',
                  '#E95C59', '#E59CC4', '#AB3282', '#23452F', '#BD956A', '#8C549C', '#585658',
                  '#9FA3A8', '#E0D4CA', '#5F3D69', '#C5DEBA', '#58A4C3', '#E4C755', '#F7F398',
                  '#AA9A59', '#E63863', '#E39A35', '#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B',
                  '#712820', '#DCC1DD', '#CCE0F5',  '#CCC9E6', '#625D9E', '#68A180', '#3A6963',
                  '#968175'
)

DimPlot(scRNA_harmony,group.by = "celltype",cols = my36colors)
save(scRNA_harmony,file = "scRNA_harmony.rdata")
########################################################################################################

scRNA1=scRNA_harmony
pt_size = 1

x=as.data.frame(scRNA1@active.ident)

library(ggplot2)
library(ggpubr)

df <- scRNA1@reductions$umap@cell.embeddings
df <- cbind(df, as.data.frame(scRNA1@active.ident))
colnames(df) <- c("x", "y", "ident")
p <- ggplot(df, aes(x=x, y=y)) + geom_point(aes(colour=ident),size=0.5) + 
  labs(x="umap_1", y ="umap_2")+ theme_pubr() +
  NoAxes() + NoLegend()
p
p1 <- p + NoLegend()
p1
DimPlot(scRNA_harmony,label = T)
p <- ggplot(df, aes(x=y, y=x)) + geom_point(aes(colour=ident),size=0.5) + 
  labs(x="umap_1", y ="umap_2")
p



p <- ggplot(df, aes(x=x, y=y)) + geom_point(aes(colour=ident),size=0.5) + 
  labs(x="umap_1", y ="umap_2")
meta = scRNA1@meta.data
meta = meta[,c(1,8)]
df.m = merge(df,meta,by=0)
colnames(df.m)
df.m <-  df.m %>%
  dplyr::group_by(seurat_clusters) %>%
  dplyr::summarise(
    x = median(x),
    y = median(y)
  )


ggplot(df, aes(x, y))  +
  geom_point(aes(colour  = ident))   +
  ggrepel::geom_text_repel(aes(label = seurat_clusters),
                           data = df.m,
                           size = 5,
                           label.size = 1,
                           segment.color = NA
  )+   theme(legend.position = "none") + theme_bw()


p+ggrepel::geom_text_repel(aes(label = seurat_clusters),
                           data = df.m,
                           size = 5,
                           label.size = 1,
                           segment.color = NA
)+   theme(legend.position = "none") + theme_bw()


### 换个主题
p+ggrepel::geom_text_repel(aes(label = seurat_clusters),
                           data = df.m,
                           size = 5,
                           label.size = 1,
                           segment.color = NA
)+   theme(legend.position = "none") + theme_dark()
##
##########################################################
p+ggrepel::geom_text_repel(aes(label = seurat_clusters),
                           data = df.m,
                           size = 5,
                           label.size = 1,
                           segment.color = NA
)+   theme(legend.position = "none") + theme_dark() + 
  labs(x="t-SNE 1", y ="t-SNE 2")+ ggtitle("huage") + ggsci::scale_color_nejm()


#########################################################################################


#########################################################################################
##自带的包也可以更改

p3 <- DimPlot(scRNA1, 
                 reduction = "umap",
                 pt.size = 0.5,
                 label = TRUE,
                 cols =  my36colors )

p3

p4 <- p3 + NoLegend()
p4

p5 = p4 + theme_dark()
p5

#################################################################################################

## 自定义绘制气泡图

B_marker= c("TCL1A","FCER2","IL4R","IGHD" , #### NaiveB
            "IFIT3","IFI44L","STAT1","ISG15",### IFN
            "HSPA1A","DNAJB1",### Activated
            "MT1X","MT2A","SLC30A1",
            "EGR1","DUSP2",#### ACB1
            "NR4A2","CD69","CD83",#### ACB2
            "CCR7","PIM3","NFKBID",
            "S100A10","CRIP1","S100A4","ITGB1","CD27","CR2","AIM2","GPR183","CD1C",
            "DUSP4","FCRL5","ZEB2","ITGAX","FGR","FCRL4","CD274",
            "NME1","PSME2","ENO1","FABP5",### PreGC
            "CXCR4","GNB2L1","ATP5L","SUGCT",### DZGC and GC
            "LMO2","GMDS","PRPSAP2","MARCKSL1",### LZGC
            "STMN1","TUBB","HMGB2","TUBA1B","MKI67",
            "JCHAIN","PRDM1","XBP1","MZB1")


library(RColorBrewer)
DotPlot(object = scRNA_harmony, features =  B_marker,scale = T,group.by = "celltype") + ##celltype_l3. ###seurat_clusters
  scale_colour_gradientn(colors=brewer.pal(9, "YlGnBu")) + theme_bw() +
  theme(axis.text.x = element_text(angle = 90)) 


B_marker1 = c("TCL1A","FCER2" ,"PRDM1","XBP1","MZB1")

f1 = sample(rownames(scRNA1),20)
p <- DotPlot(scRNA1, features = B_marker1,group.by = "celltype") + coord_flip()
p
p + RotatedAxis()


dot.data = p[["data"]]
colnames(dot.data)
p1 = ggplot(dot.data,aes(x = id, y = features.plot, color = avg.exp.scaled, size = pct.exp)) + 
  geom_point() 
p1

## 修改颜色、主题、坐标轴
ggplot(dot.data,aes(x=id, y = features.plot, color = avg.exp.scaled, size = pct.exp)) + 
  geom_point() +
  scale_color_viridis_c(name = 'expression') + 
  cowplot::theme_cowplot() + 
  theme(axis.line  = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=0.5)) +
  ylab('') +
  theme(axis.ticks = element_blank()) 

## 设置图例颜色的范围
## 使用scale_color_gradientn设置渐变色的范围为0-3，如果有的点超过了3，那么通过oob = scales::squish将它硬生生压到3
## install.packages("viridis")
p4=ggplot(dot.data,aes(x=id, y = features.plot, color = avg.exp.scaled, size = pct.exp)) + 
  geom_point() + 
  cowplot::theme_cowplot() + 
  theme(axis.line  = element_blank()) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  ylab('') +
  theme(axis.ticks = element_blank()) +
  scale_color_gradientn(colours = viridis::viridis(20), limits = c(0,3), oob = scales::squish, name = 'expression')
p4

## ggplot没有内置聚类树，但可以有以下几个选择：
# remotes::install_github("YuLab-SMU/tidytree")
#

#####################################################

## featurePlot的自定义绘制

features=c("XBP1")

VlnPlot(scRNA1, features = features ,pt.size = 0)

VlnPlot(scRNA1, features = features ,pt.size = 0, group.by = "celltype")


## 更改小提琴图横坐标顺序
Idents(scRNA1)= "celltype"

table(Idents(scRNA1))
My_levels <- c("T_cells" , "B_cell", "Monocyte","Macrophage",
                  "Fibroblasts", "Endothelial_cells","Chondrocytes", "Tissue_stem_cells", "MSC")
Idents(scRNA1) <- factor(Idents(scRNA1), levels= My_levels)
VlnPlot(scRNA1, features = features[1])

## 添加主题
FeaturePlot(scRNA1, features = features[1]) + theme_dark()
FeaturePlot(scRNA1, features = features[1])

# 加边框
FeaturePlot(scRNA1,features = features[1]) +
  annotate(geom = 'segment', y = Inf, yend = Inf, color = 'black', x = -Inf, xend = Inf, size = 1)+
  annotate(geom = 'segment', x = Inf, xend = Inf, color = 'black', y = -Inf, yend = Inf, size = 0.5) 

p = FeaturePlot(scRNA1,features = features[1],
              cols = c("lightgrey" ,"#DE1F1F"),
              slot = "data", label.size = 6, pt.size = 1.2) +  ### blank 空的
  theme(axis.line = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank())  
p

p = FeaturePlot(scRNA1,features = features[1],
              cols = c("lightgrey" ,"#DE1F1F"),
              slot = "data",label.size = 6,pt.size = 1.2) +  ### blank 空的
  theme(  axis.text = element_blank(), axis.ticks = element_blank(), axis.title = element_blank())  
p2 = p + annotate(geom = 'segment', y = Inf, yend = Inf, color = 'black', x = -Inf, xend = Inf, size = 1)+
  annotate(geom = 'segment', x = Inf, xend = Inf, color = 'black', y = -Inf, yend = Inf, size = 0.5) 
p2

## 有时候想移除图例，使用 guides()。
p + guides(fill=FALSE)
p + guides(color=F)

## 可以用theme(legend.position=…)将图例移到图表的上方、下方、左边和右边。
p2 + theme(legend.position="top")
p2 + NoLegend()

## 可以将通过rgb颜色改个透明度
print(FeaturePlot(scRNA1,features = features[1],
                  cols = c(rgb(220,212,213,180,maxColorValue = 255),
                           rgb(174,27,52,50, maxColorValue = 255)),
                  slot = "data",label.size = 6,pt.size = 1.5,
                  min.cutoff = "q10", max.cutoff = "q90") + 
        theme(axis.line = element_blank(), axis.text = element_blank(),
              axis.ticks = element_blank(), axis.title = element_blank()))



FeaturePlot(scRNA1,features = "XBP1", pt.size = 1)+
  scale_x_continuous("")+scale_y_continuous("")+
  theme_bw()+ 
  theme(  
    panel.grid.major = element_blank(),panel.grid.minor = element_blank(),  
    axis.ticks = element_blank(),axis.text = element_blank(),  
    legend.position = "none",  
    plot.title = element_text(hjust = 0.5,size=14)  
  )

#### 
exp=GetAssayData(scRNA1,slot = "data")
mat1=as.data.frame(exp["XBP1",])
colnames(mat1)="exp"
mat2=Embeddings(scRNA1,"umap")
mat3=merge(mat2,mat1,by="row.names")

colnames(mat3)
mat3%&gt;%ggplot(aes(umap_1,umap_2))+geom_point(aes(color=exp))+
  scale_color_gradient(low = "grey",high ="red")+theme_bw()+ 
  theme(  
    panel.grid.major = element_blank(),panel.grid.minor = element_blank(),  
    axis.ticks = element_blank(),axis.text = element_blank(),  
    legend.position = "none",  
    plot.title = element_text(hjust = 0.5,size=14)  
  )+
  scale_x_continuous("")+scale_y_continuous("")+labs(title= "Rbp1")


#################################################################################################
##scRNA1@active.ident=scRNA1@meta.data$celltype

## 制作marker基因的热图
my36colors <- c('#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3', '#476D87',
                  '#E95C59', '#E59CC4', '#AB3282', '#23452F', '#BD956A', '#8C549C', '#585658',
                  '#9FA3A8', '#E0D4CA', '#5F3D69', '#C5DEBA', '#58A4C3', '#E4C755', '#F7F398',
                  '#AA9A59', '#E63863', '#E39A35', '#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B',
                  '#712820', '#DCC1DD', '#CCE0F5',  '#CCC9E6', '#625D9E', '#68A180', '#3A6963',
                  '#968175'
)


Idents(scRNA1)= "celltype"
DimPlot(scRNA1)


Idents(scRNA1)= "celltype"
degs <- FindAllMarkers(scRNA1, logfc.threshold = 0.5,
                          test.use = "roc", 
                          return.thresh = 0.25, 
                          min.pct = 0.5, only.pos = T) 

degs_sig <-  degs %>% 
  filter(pct.1 > 0.3 &
           power > 0.25) %>%
  filter(cluster != "other") %>%
  arrange(cluster, -power)  

# select degs for heatmap
degs_top50 <- degs_sig %>% 
  #  filter(cluster!="other") %>%
  group_by(cluster) %>% 
  top_n(50, power) %>%
  top_n(50, avg_diff) %>%
  arrange(cluster, -power)

## 提取什么数据
exp = GetAssayData(scRNA1, slot = "data")
avgData <- exp[degs_top50$gene,] %>% 
  apply(1, function(x) {
    tapply(x, scRNA1$celltype, mean) # ExpMean
  }) %>% t


phData <- MinMax(scale(avgData), -2, 2) # z-score
rownames(phData) <- 1:nrow(phData)

##install.packages("pheatmap")
library(pheatmap)

phres <- pheatmap(
  phData, 
  color = colorRampPalette(c("darkblue", "white", "red3"))(99), #配色
  scale = "row",
  cluster_rows = F, #不按行聚类
  cluster_cols = F, #按列聚类
  clustering_method = "complete",
  show_rownames = F, #显示cluster名
  annotation_row = data.frame(cluster = degs_top50$cluster), 
)  

My_levels <- c("Macrophage" ,"Tissue_stem_cells", "Endothelial_cells",  "Monocyte",
                  "Chondrocytes","B_cell","MSC","T_cells",
                  "Fibroblasts")
phData1= phData[, My_levels ] 
phres <- pheatmap(
  phData1, 
  color = colorRampPalette(c("darkblue", "white", "red3"))(99), #配色
  scale = "row",
  cluster_rows = F, #不按行聚类
  cluster_cols = F, #按列聚类
  clustering_method = "complete",
  show_rownames = F, #显示cluster名
  annotation_row = data.frame(cluster = degs_top50$cluster), 
) 

phres
########################################################################################################
###################################################################################################

## marker基因的热图制作
### ！！！

scRNA1@active.ident=factor(scRNA1@active.ident, levels = My_levels )
table(scRNA1@active.ident)


marker = FindAllMarkers(scRNA1, logfc.threshold = 0.5, min.pct = 0.5, only.pos = T)

write.csv(marker,file = "marker.csv")

top10 <-  marker  %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

scRNA1 <-  ScaleData(scRNA1, features = rownames(scRNA1))
DoHeatmap(scRNA1,features = top10$gene, label = F, slot = "scale.data")
table(scRNA1@active.ident)

sc.s = subset(scRNA1,downsample=500)
DoHeatmap( sc.s ,features = top10$gene,label = F,slot = "scale.data")
?DoHeatmap

table(scRNA1$celltype)
scRNA2=scRNA1

table(scRNA1@active.ident)
scRNA1 = subset(scRNA1, downsample = 150)
library(pheatmap)

colanno = scRNA1@meta.data 
colanno$barcode = rownames(colanno)
colanno = colanno%>%arrange(celltype)
rownames(colanno) = colanno$barcode
colanno$barcode = NULL



colanno1 = colanno[,10]
colanno1 = as.data.frame(colanno1)
rownames(colanno1) = rownames(colanno)
colnames(colanno1) = "celltype"

rowanno=top10
colnames(top10)
rowanno=rowanno%>%arrange(cluster)

scale.data=GetAssayData(scRNA1,slot =  "scale.data")
mat4=scale.data[rowanno$gene,rownames(colanno1)]
mat4[mat4 >= 2.5]=2.5
mat4[mat4 < (-1.5)]= -1  

pheatmap(mat4,cluster_rows = F,cluster_cols = F,
         show_colnames = F,
         annotation_col = colanno1,
         gaps_row=as.numeric(cumsum(table(rowanno$cluster))[-9]),
         gaps_col=as.numeric(cumsum(table(colanno$celltype))[-9]) 
         
)



pheatmap(mat4,cluster_rows = F,cluster_cols = F,
         show_colnames = F,
         annotation_col = colanno1,
         gaps_row=as.numeric(cumsum(table(rowanno$cluster))[-6]),
         gaps_col=as.numeric(cumsum(table(colanno$celltype))[-6]),
         color = colorRampPalette(c("darkblue", "white", "red3"))(99)
         
)



pheatmap(mat4,cluster_rows = F,cluster_cols = F,
         show_colnames = F,
         annotation_col = colanno1,
         gaps_row=as.numeric(cumsum(table(rowanno$cluster))[-6]),
         gaps_col=as.numeric(cumsum(table(colanno$celltype))[-6]) ,
         filename="marker.heatmap.pdf",width=11,height = 7
)


################################################################################################


######################################################################################################

## 自定义绘制一个小提琴图

top3 <-  marker %>% group_by(cluster) %>% top_n(n = 3, wt = avg_log2FC)

scRNA1 = scRNA_harmony
VlnPlot(scRNA1, features = top3$gene[1:3] , pt.size = 0, ncol = 1)+
  scale_x_discrete("")+
  theme(
    axis.text.x.bottom = element_blank()
  )


Idents(scRNA1) = "celltype"
marker = FindAllMarkers(scRNA1,logfc.threshold = 0.5, min.pct = 0.5, only.pos = T)

top10 <- marker %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)


library(reshape2)
exp = GetAssayData(scRNA1,slot = "data")
vln.df = as.data.frame(exp[top10$gene,])
vln.df$gene = rownames(vln.df)
vln.df = melt(vln.df,id="gene")
colnames(vln.df)[c(2,3)] = c("barcode","exp")

anno = scRNA1@meta.data
anno$barcode = rownames(anno)


vln.df = inner_join(vln.df,anno,by="barcode")
##vln.df$gene=factor(vln.df$gene,levels = top10$gene[1:7]) #为了控制画图的基因顺序
table(vln.df$gene)
vln.df$gene = factor(vln.df$gene,levels = c(.........))

vln.df%>%ggplot(aes(celltype,exp))+geom_violin(aes(fill=gene),scale = "width")+
  facet_grid(vln.df$gene~.,scales = "free_y")+
  scale_fill_brewer(palette = "Set3",direction = 1)+
  scale_x_discrete("")+scale_y_continuous("")+
  theme_classic()+
  theme(
    axis.text.x.bottom = element_text(angle = 45,hjust = 1,vjust = 1),
    panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
    legend.position = "none"
  )



###########################################################################################################
## 


## 
top10 <- marker  %&>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
DotPlot(scRNA1, features = top10$gene,cols = c("lightgrey", "#FF00FF"))+
  scale_x_discrete("huage")+scale_y_discrete("")+theme_classic()+RotatedAxis()
# 
?DotPlot
bubble.df = as.matrix(scRNA1[["RNA"]]@data[ top10$gene,])
bubble.df = t(bubble.df)
bubble.df = as.data.frame(scale(bubble.df))
bubble.df$CB = rownames(bubble.df)
meta = scRNA1@meta.data[,c(1,6)]
bubble.df = merge(bubble.df,meta,by.x = "CB",by.y = 0)
bubble.df$CB = NULL

celltype_v=c()
gene_v=c()
mean_v=c()
ratio_v=c()
for (i in unique(bubble.df$celltype)) {
  bubble.df_small=bubble.df%&gt;%filter(celltype==i)
  for (j in top10$gene) {
    exp_mean=mean(bubble.df_small[,j])
    exp_ratio=sum(bubble.df_small[,j] &gt; min(bubble.df_small[,j])) / length(bubble.df_small[,j])
    celltype_v=append(celltype_v,i)  ##Add elements to a vector.
    gene_v=append(gene_v,j)
    mean_v=append(mean_v,exp_mean)
    ratio_v=append(ratio_v,exp_ratio)
  }
}

plotdf=data.frame(
  celltype=celltype_v,
  gene=gene_v,
  exp=mean_v,
  ratio=ratio_v
)

library(RColorBrewer)
plotdf$celltype=factor(plotdf$celltype,levels = sort(unique(plotdf$celltype)))
plotdf$gene=factor(plotdf$gene,levels = rev(as.character(top10$gene[1:10])))
plotdf$exp=ifelse(plotdf$exp&gt;3,3,plotdf$exp)
plotdf%>%ggplot(aes(x=celltype,y=gene,size=ratio,color=exp))+geom_point()+
  scale_x_discrete(" ")+scale_y_discrete(" ")+
  scale_color_gradientn(colours = rev(c("#FFD92F","#FEE391",brewer.pal(11, "Spectral")[7:11])))+
  scale_size_continuous(limits = c(0,1))+theme_classic()+
  theme(
    axis.text.x.bottom = element_text(hjust =1, vjust = 1, angle = 90)
  )






################################################################################


library(ggpubr)

VlnPlot( scRNA1 , "XBP1",group.by = "orig.ident")+ 
  theme_bw()+stat_compare_means(method = 't.test')

VlnPlot( scRNA1 , "XBP1",group.by = "orig.ident")+ 
  theme_bw()+stat_compare_means(method = 't.test')+geom_boxplot()

VlnPlot( scRNA1 , "XBP1",group.by = "celltype")+ 
  theme_bw()+stat_compare_means(method = 't.test')+geom_boxplot()

VlnPlot(subset(scRNA1,XBP1 &gt; 0 ), "XBP1",group.by = "orig.ident")+ 
  theme_bw()+stat_compare_means(method = 't.test')+geom_boxplot()

sc.FURIN=subset(scRNA1,Rbp1 &gt; 0 )
FURIN.exp=sc.FURIN@assays$RNA@data["Rbp1",]
FURIN.exp.m= sc.FURIN@meta.data
table(sc.FURIN$orig.ident)
FURIN.exp.m$exp=FURIN.exp
FURIN.exp.m$Cohort=factor(FURIN.exp.m$orig.ident,
                          levels = c("DapiNeg1","DapiNeg2","CAF1",'CAF2'))
write.csv(FURIN.exp.m,file = "FURIN.exp.m.csv")
colnames(FURIN.exp.m)
library(ggplot2)
library("ggsignif")
my_comparisons=list(c("DapiNeg1","CAF1"),c("DapiNeg2",'CAF2'),
                    c("DapiNeg2" ,"CAF1"))

ggplot(FURIN.exp.m,aes(Cohort,exp,fill=Cohort)) +
  geom_violin()+theme_classic()+
  theme(plot.title=element_text(size = 15,face="bold"),
        axis.text.x=element_text(size=15,face="bold",angle=25,hjust = 1,vjust = 1),
        axis.text.y=element_text(size=18,face="bold"),
        axis.title.x=element_text(size = 0,vjust = 1,face="bold"),
        axis.title.y=element_text(size = 15,face="bold"))+
  labs(y= 'Expression',x="Group")+
  geom_signif(comparisons = my_comparisons,step_increase = 0.1,
              map_signif_level = F,test = t.test,size=1,textsize = 6)+
  geom_jitter(shape=16, position=position_jitter(0.2))+
  ggtitle("Rbp1")+
  theme(plot.title = element_text(size =18,hjust = 0.5))+
  guides(fill=guide_legend(title = " "))+ ##更改legend名字
  ##如果是color分组 需要 guides(color=guide_legend(title = "Group")) 
  theme(axis.line.x=element_line(linetype=1,color="black",size=1))+
  theme(axis.line.y=element_line(linetype=1,color="black",size=1))+
  theme(legend.text = element_text(size = 15) ) ##legend字体变大

差异基因分析

差异基因的展示（heatmap等）

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
## 差异基因分析

##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())

###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(harmony)
library(cowplot)
library(ggrepel)
library(reshape2)

setwd("D:/0_SubjectData/")
load("scRNA_harmony.rdata")

table(scRNA_harmony@active.ident)
setwd("D:/0_SubjectData/")

#######################################################################################
library(ggplot2)
library(cowplot)
theme_set(theme_cowplot())

t.cells <- subset(scRNA_harmony, idents = "1")

table(t.cells$orig.ident)

## 对两个分组进行差异分析
## ident.1，ident.2
deg <- FindMarkers(t.cells, ident.1 = "sample_11",ident.2 = "sample_3", 
                      logfc.threshold =1.5, group.by = "orig.ident")


Idents(t.cells) = "orig.ident"
deg <- FindMarkers(t.cells, ident.1 = "sample_11", ident.2 = "sample_3", 
                      logfc.threshold = 1.5)



exp = AverageExpression(t.cells, verbose = FALSE)
exp = as.data.frame(exp$RNA)  
## 取log值是防止基因的离散程度过高，表达量太大
exp = log1p(exp)


##avg.t.cells <- as.data.frame(log1p(AverageExpression(t.cells, verbose = FALSE)$RNA))
avg.t.cells = exp

avg.t.cells$gene <- rownames(avg.t.cells)
colnames(avg.t.cells)

deg <- FindMarkers(t.cells,ident.1 = "sample_11",ident.2 = "sample_3", 
                      logfc.threshold =8, group.by = "orig.ident")


deg <- FindMarkers(t.cells,ident.1 =c("sample_11","sample_12") ,
                      ident.2 = c("sample_3","sample_4") , 
                      logfc.threshold =8, group.by = "orig.ident")

colnames(avg.t.cells)
colnames(avg.t.cells) = c("sample_11","sample_3","gene")

genes.to.label = sample(rownames(deg),20)
p1 <- ggplot(avg.t.cells, aes(sample_11,sample_3)) + geom_point() + ggtitle("T Cells")
p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p1 

p1 <- ggplot(avg.t.cells, aes(sample_11,sample_3)) + geom_point() + ggtitle("T Cells")
LabelPoints(plot = p1, points = c("RPS4Y1"), repel = TRUE)


###########################################################

degdf = avg.t.cells
degdf.l = degdf[genes.to.label,]
degdf.p = degdf[-match(genes.to.label,rownames(degdf)),]
## 标注为不显著
degdf.p$Sig = "NO_DIFF"
## 标注为显著
degdf.l$Sig = "DIFF"
degdf = rbind(degdf.p,degdf.l)
colnames(degdf)
p = ggplot(degdf, aes(sample_11,sample_3)) +
  geom_point( aes(color=Sig)) +
  
  scale_color_manual(values=c("red","grey"))+
  ggtitle("Stromal Cells")+
  theme(plot.title = element_text(size =18,hjust = 0.5, face = "bold")) +
  
  theme(axis.title.x =element_text(size=14), axis.title.y=element_text(size=14)) + 
  
  theme_bw()

p

p = p+theme_bw()+
  theme(panel.border = element_blank(),panel.grid.major = element_blank(),
        
        panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"))

p
data_selected <- degdf[genes.to.label,]
p + geom_label_repel(data=data_selected,
                     aes(label=rownames(data_selected)))


?geom_point
p +
  geom_point(size = 6, shape =2, data = data_selected[c("MMP2"),],colour="green") +
  ggrepel::geom_text_repel(
    aes(label =genes.to.label),size =2.5,
    data = data_selected,
    color="blue")+theme_bw()+
  theme(panel.border = element_blank(),panel.grid.major = element_blank(),
        
        panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"))+
  ggtitle("Stromal Cells")+
  theme(plot.title = element_text(size =18,hjust = 0.5))


?geom_label_repel


#######################################################################################
##火山图
table(scRNA_harmony@active.ident)
cd4.naive=subset(scRNA_harmony,idents ="1")

degdf <- FindMarkers(cd4.naive,ident.1 = "sample_11",ident.2 = "sample_3", logfc.threshold = 0.01,group.by = "orig.ident")

degdf$symbol <- rownames(degdf)
logFC_t = 0
P.Value_t = 1e-28
degdf$change = ifelse(degdf$p_val_adj < P.Value_t & degdf$avg_log2FC < 0,"down",
                      ifelse(degdf$p_val_adj < P.Value_t & degdf$avg_log2FC > 0,"up","stable"))
p=ggplot(degdf, aes(x = avg_log2FC, y = -log10(p_val_adj), size = pct.1)) +
  geom_point(alpha = 0.4,  aes(color = change)) +
  ylab("-log10(Pvalue)")+
  scale_color_manual(values=c("green", "grey", "red"))+
  geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
  theme_bw()


## 
p=p+theme_bw()+
  theme(panel.border = element_blank(),panel.grid.major = element_blank(),
        
        panel.grid.minor = element_blank(),axis.line = element_line(colour = "black"))


## 某个基因的单独展示
data_selected <- degdf["TIMP1",]
p + geom_label_repel(data=data_selected,
                     aes(label=rownames(data_selected)))


p + geom_text_repel(data=data_selected,
                    aes(label=rownames(data_selected)))

?geom_point
p +
  geom_point(size = 4, shape =2, data = data_selected,colour="blue") +
  ggrepel::geom_label_repel(
    aes(label = symbol),
    data = data_selected,
    color="#f7aa5d")
###如何改变点的颜色

#############################################################################################

三维pca的展示

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163

## 三维pca的展示

##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())

library(Seurat)
library(scatterplot3d)

setwd("D:/0_SubjectData/")
load("sc.combined.rdata")
setwd("D:/0_SubjectData/")

table(sc.combined@meta.data$new.celltype.2)
table(sc.combined@meta.data$type)

table(sc.combined@meta.data$orig.ident)


### 三维pca图可以直接观看不同亚群之间的差异性
##  Fibroblasts

sc.imm=subset(sc.combined,idents = c("Epithelial","Endothelial","Stem_cell"))


Idents(sc.imm)="orig.ident"
table( sc.imm@active.ident) 
# factor 可以进行排序
sc.imm@meta.data$orig.ident=factor(sc.imm@meta.data$orig.ident,levels = c(
  "GF_1","GF_2","GF_3","SPF_1","SPF_2","SPF_3","FMT_1","FMT_2","FMT_3"
))
Idents(sc.imm)="orig.ident"
table(sc.imm@active.ident)
?AverageExpression
x=AverageExpression(sc.imm ,add.ident="new.celltype.2")
x1=AverageExpression(sc.imm  )
x=x$RNA
x=as.matrix(x)
x1=x1$RNA
x1=as.matrix(x1)
pca <- prcomp(t(x))

pca.data=pca$x

rownames(pca.data)

pca.data=as.data.frame(pca.data)
pca.data=pca.data[,1:3]

pca.data$Type = c(rep("GF",9),rep("SPF",9),rep("FMT",9))

s1=strsplit(rownames(pca.data),split = "_",fixed = T)
type=sapply(s1, function(x){x[1]}   )


pca.data$cell.type = c(rep(c("Epithelial","Endothelial","Stem_cell"),9 ))



Type.p=c(rep(11,9),rep(16,9),rep(17,9))
cell.type.p = c(rep(c("blue","red","orange"),9 ))


colors.lib <- c("blue","red","orange")
shapes.lib = c(11,16,17)

getwd()

# 1. Source the function
source( "E:/super.lesson/lesson9/huage.3D.plot.R" )
# 2. 3D scatter plot


s3d <- scatterplot3d(pca.data[,c("PC1","PC2","PC3")],
                        
                        pch = Type.p,color = cell.type.p,
                        
                        cex.symbols = 1,grid=FALSE, box=FALSE, 
                        
                        main = "3D PCA plot")


legend("topright",
       c("Epithelial","Endothelial","Stem_cell"),
       fill=c('blue',"red","orange"),
       box.col=NA)

legend("topleft",
       c('GF','SPF','FMT'),
       pch = shapes.lib,
       box.col=NA)

# 3. Add grids
addgrids3d(pca.data[,c("PC1","PC2","PC3")], grid = c("xy", "xz", "yz"))

table(sc.combined@meta.data$new.celltype.2)

##  Goblet Paneth cells Enteroendocrine
Idents(sc.combined)="new.celltype.3"
sc.imm=subset(sc.combined,idents = c("Goblet","MSC","T_NK"))


Idents(sc.imm)="orig.ident"
x=AverageExpression(sc.imm ,add.ident="new.celltype.3")
x=x$RNA
pca <-  prcomp(t(x))

pca.data=pca$x

rownames(pca.data)

pca.data=as.data.frame(pca.data)
pca.data=pca.data[,1:3]

pca.data$Type = c(rep("GF",9),rep("SPF",9),rep("FMT",9))
pca.data$cell.type = c(rep(c("Goblet","Paneth cells","Enteroendocrine"),9 ))



Type.p=c(rep(11,9),rep(16,9),rep(17,9))
cell.type.p = c(rep(c("blue","red","orange"),9 ))


colors.lib &lt;- c("blue","red","orange")
shapes.lib = c(11,16,17)


getwd()


# 1. Source the function
source('plot.3d.R')
# 2. 3D scatter plot


s3d &lt;- scatterplot3d(pca.data[,c("PC1","PC2","PC3")],
                        
                        pch = Type.p,color = cell.type.p,
                        
                        cex.symbols = 1,grid=FALSE, box=FALSE, 
                        
                        main = "3D PCA plot")




legend("topright",
       c("Goblet","Paneth cells","Enteroendocrine"),
       fill=c('blue',"red","orange"),
       box.col=NA)

legend("topleft",
       c('GF','SPF','FMT'),
       pch = shapes.lib,
       box.col=NA)

# 3. Add grids
addgrids3d(pca.data[,c("PC1","PC2","PC3")], grid = c("xy", "xz", "yz"))

多个细胞类型基因展示

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
## 不同分组差异基因的展示

##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())
library(dplyr)
library(Seurat)
library(patchwork)
library(reshape2)
library(RColorBrewer)
library(ggplot2)
library(ggrepel)  
library(magrittr)
library(data.table)
library(dplyr)
library(Seurat)
library(tidyverse)
library(patchwork)

setwd("D://0_SubjectData/")
load("scRNA_harmony.rdata")

DimPlot(scRNA_harmony,group.by = "orig.ident")

DimPlot(scRNA_harmony,group.by = "celltype")


setwd("D://0_SubjectData/")
sc.combined=scRNA_harmony

table(sc.combined@meta.data$celltype)
table(sc.combined@meta.data$orig.ident)
type=c("Chondrocytes","Endothelial_cells","Fibroblasts",
       "Macrophage","Monocyte", "T_cells", 
       "Tissue_stem_cells")

## 等价于上面的代码
type = unique(sc.combined@meta.data$celltype)

dir.create("deg/")
setwd("deg/")




#library(future)
#plan("multiprocess", workers =5)
#options(future.globals.maxSize = 2000 * 1024^2)

r.deg = data.frame()
table(sc.combined@meta.data$orig.ident)

for (i in 1:length(type)) {
  Idents(sc.combined)="celltype"
  deg = FindMarkers(sc.combined,ident.1 = "sample2",ident.2 = "sample21",
                  group.by = "orig.ident",subset.ident =type[i]   )
  
  write.csv(deg,file = paste0( type[i],'deg.csv') )
  deg$gene = rownames(deg)
  deg$celltype = type[i]
  deg$unm = i-1
  r.deg = rbind(deg, r.deg)
  
}

Idents(sc.combined)="celltype"
deg=FindMarkers(sc.combined,ident.1 = "sample_11",ident.2 = "sample_3",
                group.by = "orig.ident",subset.ident =type[2]   )


table(r.deg$unm) 
#############################################################################

# 根据自己计算的marker基因数量确定log2FC的阈值，这里先定为1.5
r.deg <- subset(r.deg, p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)
r.deg$threshold <- as.factor(ifelse(r.deg$avg_log2FC > 0 , 'Up', 'Down'))
dim(r.deg)

r.deg$adj_p_signi <- as.factor(ifelse(r.deg$p_val_adj < 0.01 , 'Highly', 'Lowly'))
r.deg$thr_signi <- paste0(r.deg$threshold, "_", r.deg$adj_p_signi)
r.deg$unm <- &>% as.vector(.) %&>% as.numeric(.)

##自定义显示想要展示的基因名
##这里挑选log2FC为top5的基因进行展示

top_up_label <- r.deg %>% 
  subset(., threshold%in%"Up") %>% 
  group_by(unm) %>% 
  top_n(n = 5, wt = avg_log2FC) %>% 
  as.data.frame()

top_down_label <- r.deg %>% 
  subset(., threshold %in% "Down") %>% 
  group_by(unm) %>% 
  top_n(n = -5, wt = avg_log2FC) %>% 
  as.data.frame()

top_label <- rbind(top_up_label,top_down_label)
top_label$thr_signi %<>% 
  factor(., levels = c("Up_Highly","Down_Highly","Up_Lowly","Down_Lowly"))

# 保存到文件，便于小白套用格式
#write.csv(top_label, "easy_input_label.csv", quote = F)
##也可以基于output_pbmc.markers.csv文件，手动挑选出想要标注名字的基因，
#例如标注参与某一通路的基因，然后将文件命名为easy_input_label.csv


colnames(r.deg)
### 准备绘制暗灰色背景所需数据 
background_position <- r.deg %>%
  dplyr::group_by(unm) %>%
  dplyr::summarise(Min = min(avg_log2FC) - 0.2, Max = max(avg_log2FC) + 0.2) %>%
  as.data.frame()
## summarise()` ungrouping output (override with `.groups` argument)
background_position$unm %<>% as.vector(.) %>% as.numeric(.)
background_position$start <- background_position$unm - 0.4
background_position$end <- background_position$unm + 0.4

### 准备绘制中间区域cluster彩色bar所需数据
cluster_bar_position <- background_position
cluster_bar_position$start <- cluster_bar_position$unm - 0.5
cluster_bar_position$end <- cluster_bar_position$unm + 0.5
cluster_bar_position$unm %<>% 
  factor(., levels = c(0:max(as.vector(.))))

## 设置填充颜色
cols_thr_signi <- c("Up_Highly" = "#d7301f",
                       "Down_Highly" = "#225ea8",
                       "Up_Lowly" = "black",
                       "Down_Lowly" = "black")
cols_cluster <- c("0" = "#35978f",
                     "1" = "#8dd3c7",
                     "2" = "#ffffb3",
                     "3" = "#bebada",
                     "4" = "#fb8072",
                     "5" = "#80b1d3",
                     "6" = "#fdb462","7" = "#925bea","8" = "#db5e92" 
)

p= ggplot() +
  geom_rect(data = background_position, aes(xmin = start, xmax = end, ymin = Min,
                                            ymax = Max),
            fill = "#525252", alpha = 0.1) + ###添加灰色背景色
  geom_jitter(data = r.deg, aes(x =unm, y = avg_log2FC, colour = thr_signi),
              size = 1,position = position_jitter(seed = 1)) +
  scale_color_manual(values = cols_thr_signi) +
  scale_x_continuous(limits = c(-0.5, max(r.deg$unm) + 0.5),
                     breaks = seq(0, max(r.deg$unm), 1),
                     label = seq(0, max(r.deg$unm),1)) + #修改坐标轴显示刻度
  
  # 根据top_label标注基因名
  geom_text_repel(data = top_label, aes(x =unm, y = avg_log2FC, label = gene),
                  position = position_jitter(seed = 1), show.legend = F, size = 2.5,
                  box.padding = unit(0, "lines")) +
  
  geom_rect(data = cluster_bar_position, aes(xmin = start, xmax = end, ymin = -0.4,
                                             ymax = 0.4, fill = unm), color = "black", alpha = 1, show.legend = F) +
  scale_fill_manual(values = cols_cluster) +
  labs(x = "Cluster", y = "average log2FC") +
  theme_bw()

plot1 <-  p + theme(panel.grid.minor = element_blank(), ##去除网格线
                      panel.grid.major = element_blank(),
                      axis.text.y = element_text(colour = 'black', size = 14),
                      axis.text.x = element_text(colour =c("#89288F","#89288F","#F47D2B","#F47D2B","#FF00FF",'black','black'), size = 14, vjust = 59.5), #调整x轴坐标,vjust的值按照最终结果稍加调整
                      panel.border = element_blank(), ## 去掉坐标轴
                      axis.ticks.x = element_blank(), ## 去掉的坐标刻度线
                      axis.line.y = element_line(colour = "black")) #添加y轴坐标轴

plot1
?element_text
ggsave(filename = "deg_pointplot.pdf", plot = plot1, width = 9, height = 6)


通路富集分析

image-20240921164624949

GO 指的是 Gene Ontology（基因本体论），是对基因在不同维度和不同层次上的描述。

GO 分为三个主要类别：

  1. Biological Process（生物过程）
    • 描述基因产物参与的生物学过程或路径，例如：细胞分裂、DNA 复制、凋亡（细胞死亡）等。
  2. Molecular Function（分子功能）
    • 描述基因产物的基本功能，通常指特定的分子活动，比如：催化活性、结合活性（如 ATP 结合）等。
  3. Cellular Component（细胞成分）
    • 描述基因产物所在的细胞部位或细胞结构，比如：线粒体、细胞核、膜等。
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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
## 通路富集分析

## 1. 处理组和对照组进行差异富集分析

## 2. 针对于相应的亚群的marker基因进行通路的富集分析


## 系统报错改为英文
Sys.setenv(LANGUAGE = "en")
## 禁止转化为因子
options(stringsAsFactors = FALSE)
## 清空环境
rm(list=ls())

setwd("/home/oelab/RunData2/heisenberg/")

### 加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(harmony)
library(cowplot)
library(ggplot2)

scRNA_harmony <- readRDS("/home/oelab/RunData2/heisenberg/scRNA_harmony.rds")

scRNA_harmony <- JoinLayers(scRNA_harmony)

## BiocManager::install("AUCell")
library(AUCell)
library(ggplot2)
library(Seurat)
## BiocManager::install("clusterProfiler")
library(clusterProfiler)


## 对细胞进行抽样
sc.id = sample(colnames(scRNA_harmony), 1500)
sc2 = scRNA_harmony[,sc.id]

## install.packages("doParallel")
## install.packages("doRNG")

## 对表达矩阵进行提取
names(sc2@assays$RNA@layers)
colnames(sc2@assays$RNA@layers$data)

# 如果基因名称缺失，可以从 Seurat 对象中提取基因名并设置为行名
rownames(sc2@assays$RNA@layers$data) <- rownames(sc2)

# 如果细胞名称缺失，可以从 Seurat 对象中提取细胞名称并设置为列名
colnames(sc2@assays$RNA@layers$data) <- colnames(sc2)

cells_rankings <- AUCell_buildRankings(sc2@assays$RNA@layers$data,  nCores=6, plotStats=TRUE) 

cells_rankings

## 读取KEGG基因集文件，是在GSEA网站下载
c2 <- read.gmt("/home/oelab/RunData2/heisenberg/AllData/c2.cp.kegg_medicus.v2024.1.Hs.symbols.gmt") 

## lapply是列表的循环，是为了得到每个基因集的列表形式
geneSets <- lapply(unique(c2$term), function(x){print(x);c2$gene[c2$term == x]})
names(geneSets) <- unique(c2$term)

## 
?AUCell_calcAUC
cells_AUC <-  AUCell_calcAUC(geneSets, cells_rankings,nCores =1,
                               aucMaxRank=nrow(cells_rankings)*0.1)


grep("OX", rownames(cells_AUC@assays@data$AUC),value = T)

geneSet <- "KEGG_MEDICUS_REFERENCE_BETA_OXIDATION" 

## 提取通路的数据
aucs <- as.numeric(getAUC(cells_AUC)[geneSet, ])
sc2$AUC <- aucs
df <-  data.frame(sc2@meta.data, sc2@reductions$umap@cell.embeddings)
colnames(df)
class_avg <- df %>%
  group_by(celltype) %>%
  summarise(
    umap_1 = median(umap_1),
    umap_2 = median(umap_2)
  )

ggplot(df, aes(umap_1, umap_2))  +
  geom_point(aes(colour  = AUC)) + viridis::scale_color_viridis(option="D") +
  ggrepel::geom_label_repel(aes(label = celltype),
                            data = class_avg,
                            size = 3,
                            label.size = 1,
                            segment.color = NA
  )+   theme(legend.position = "none") + theme_bw() + facet_grid(.~orig.ident)

colnames(df)

library(viridis)
ggplot(data.frame(sc2@meta.data,sc2@reductions$umap@cell.embeddings), aes(umap_1, umap_2, color=AUC)
) + geom_point( size=1.5
) + scale_color_viridis(option="A")  + theme_light(base_size = 26) + facet_grid(.~orig.ident)


## 文章中的代码
ggplot(data.frame(cd8.seurat@meta.data, cd8.seurat@dr$umap@cell.embeddings), aes(UMAP1, UMAP2, color=AUC)
) + geom_point( size=1.5
) + scale_color_viridis(option="A")  + theme_light(base_size = 26) + facet_grid(.~Type)


## 直接调用msigDB的数据库
## install.packages("msigdbr")
library(msigdbr)
msigdbr_species()
msigdbr_collections()

m_df <- msigdbr(species = "Homo sapiens",  category = "C2", subcategory = "KEGG")
fgsea_sets <-  m_df %>% split(x = .$gene_symbol, f = .$gs_name)



#############################################################################################
###  gsva

### 通路的打分算法 GSVA
### 是基于基因集的分析方法，主要用于表征不同样本中基因表达的变化，而不是单个基因的表达水平

library('GSEABase')
library(GSVA)
library(msigdbr)
m_df <- msigdbr(species ="Homo sapiens",  category = "H" )
geneSets <- m_df %>% split(x = .$gene_symbol, f = .$gs_name)
colnames(scRNA_harmony@meta.data)
table(scRNA_harmony$celltype)
Idents(scRNA_harmony)="celltype"
## 把各个类型的细胞取平均值
?AverageExpression
exp=AverageExpression(scRNA_harmony) 
exp=exp[["RNA"]]
exp=as.matrix(exp)

?gsva
?gsvaParam

# 创建 gsvaParam 对象
gsva_params <- gsvaParam(
  exprData = as.matrix(exp),  # 输入表达矩阵
  geneSets = geneSets,         # 基因集
  minSize = 5,                 # 最小基因集大小
  maxSize = 500,               # 最大基因集大小
  kcdf = "Gaussian",           # 使用高斯分布
  tau = 1,                     # tau 参数
  maxDiff = TRUE,              # 使用最大差异排名
  absRanking = FALSE            # 不使用绝对排名
)

# 调用 gsva 函数
GSVA_hall <- gsva(
  param = gsva_params,
  verbose = TRUE,
)

# GSVA_hall <- gsva(expr=as.matrix(exp), 
#                      gset.idx.list=geneSets, 
#                      mx.diff=T, # 数据为正态分布则T，双峰则F
#                      kcdf="Gaussian", #CPM, RPKM, TPM数据就用默认值"Gaussian"， read count数据则为"Poisson"，
#                      parallel.sz=4) # 并行线程数目



head(GSVA_hall)
?pheatmap::pheatmap
p <- pheatmap::pheatmap(GSVA_hall, #热图的数据
                   # GSVA_hall[1:20, ],
                   cluster_rows = T,#行聚类
                   cluster_cols =T,#列聚类，可以看出样本之间的区分度
                   
                   show_colnames=T,
                   scale = "row", #以行来标准化，这个功能很不错
                   color =colorRampPalette(c("blue", "white", "red"))(100))

ggsave("gsva_pheatmap.pdf", p, width = 8,height  = 8)

##################################################################################


## 对不同的分组对象进行GSVA打分

exp=AverageExpression(scRNA_harmony,add.ident = "orig.ident") 
exp=exp[["RNA"]]

counts2=exp[,c(1:6)]

# 创建 gsvaParam 对象
gsva_params <- gsvaParam(
  exprData = as.matrix(exp),  # 输入表达矩阵
  geneSets = geneSets,         # 基因集
  minSize = 5,                 # 最小基因集大小
  maxSize = 500,               # 最大基因集大小
  kcdf = "Gaussian",           # 使用高斯分布
  tau = 1,                     # tau 参数
  maxDiff = TRUE,              # 使用最大差异排名
  absRanking = FALSE            # 不使用绝对排名
)

# 调用 gsva 函数
GSVA_hall <- gsva(
  param = gsva_params,
  verbose = TRUE,
)

head(GSVA_hall)

pheatmap::pheatmap(GSVA_hall, #热图的数据
                   cluster_rows = T,#行聚类
                   cluster_cols =F,#列聚类，可以看出样本之间的区分度
                   
                   show_colnames=T,
                   scale = "row", #以行来标准化，这个功能很不错
                   color =colorRampPalette(c("#FF7744", "white","#AAAAAA","#0044BB"))(100))

#################################################################################################
#################################################################################################
############################################################################################


### 绘制CNS上面的通路分析图示

##install.packages("msigdbr")


##BiocManager::install("GSVA")
library('GSEABase')
library(GSVA)
library(msigdbr)
m_df <- msigdbr(species ="Homo sapiens",  category = "H" )
geneSets <- m_df %>% split(x = .$gene_symbol, f = .$gs_name)

table(scRNA_harmony$celltype ,scRNA_harmony$orig.ident)

Idents(scRNA_harmony)="celltype"
sc.b = subset(scRNA_harmony,ident="Micro.")
table(sc.b$orig.ident)
Idents(sc.b)="orig.ident"

# 定义处理组和对照组的标识（这里定义了一个包含“处理组样本标识的向量”和一个包含“对照组样本标识的向量”。）  
treatments <- c("Human_Organoid_AD","Human_3D_AD", "Human_AD_01","Human_AD_02","Human_AD_03","Human_AD_04","Human_AD_05","Human_AD_06",
                "Human_AD_07","Human_AD_08","Human_AD_09","Human_AD_10")  
controls <- c("Human_Organoid_Normal","Human_3D_Normal", "Human_Normal_01","Human_Normal_02","Human_Normal_03","Human_Normal_04","Human_Normal_05",
              "Human_Normal_06","Human_Normal_07")   

# 创建一个新的标识列，将处理组合并为一个组，对照组合并为一个组  
sc.b[["treat_vs_control"]] <- 
  ifelse(sc.b$orig.ident %in% treatments, #判断每个样本的原始标识（orig.ident）是否在处理组的标识中
         "Treated",
         "Control")  

# 3. 将"treat_vs_control"设置为astrocyte_subset新的标识
Idents(sc.b) <- "treat_vs_control"

sc.b = subset(sc.b, ident=c("Treated", "Control"))

sc.b <- JoinLayers(sc.b)
table(sc.b$orig.ident)
exp=GetAssayData(sc.b,slot = "data")
exp=as.matrix(exp)


# 创建 gsvaParam 对象
gsva_params <- gsvaParam(
  exprData = as.matrix(exp),  # 输入表达矩阵
  geneSets = geneSets,         # 基因集
  minSize = 5,                 # 最小基因集大小
  maxSize = 500,               # 最大基因集大小
  kcdf = "Gaussian",           # 使用高斯分布
  tau = 1,                     # tau 参数
  maxDiff = TRUE,              # 使用最大差异排名
  absRanking = FALSE            # 不使用绝对排名
)

# 调用 gsva 函数
GSVA_hall <- gsva(
  param = gsva_params,
  verbose = TRUE,
)

head(GSVA_hall)

## limma差异通路分析

# BiocManager::install('limma')
library(limma)
# 设置或导入分组

group <- factor(sc.b@meta.data$orig.ident, levels = c( 'sample2','sample21'))
design <- model.matrix(~0+group)
colnames(design) = levels(factor(group))
rownames(design) = colnames(GSVA_hall)
design
# Tunor VS Normal
compare <- makeContrasts(sample2 - sample21, levels=design)
fit <- lmFit(GSVA_hall, design)
fit2 <- contrasts.fit(fit, compare)
fit3 <- eBayes(fit2)
Diff <- topTable(fit3, coef=1, number=200)
head(Diff)

## 发散条形图绘制
## barplot
dat_plot <- data.frame(id = row.names(Diff),
                          t = Diff$t)
# 去掉"HALLMARK_"
library(stringr)
dat_plot$id <- str_replace(dat_plot$id , "HALLMARK_","")
# 新增一列 根据t阈值分类
dat_plot$threshold = factor(ifelse(dat_plot$t  > -1, ifelse(dat_plot$t >= 1 ,'Up','NoSignifi'),'Down'),levels=c('Up','Down','NoSignifi'))
# 排序
dat_plot <-  dat_plot %&>% arrange(t)
# 变成因子类型
dat_plot$id <- factor(dat_plot$id,levels = dat_plot$id)
# 绘制
library(ggplot2)

##install.packages("ggthemes")
library(ggthemes)
# install.packages("ggprism")
library(ggprism)
p <- ggplot(data = dat_plot,aes(x = id,y = t,fill = threshold)) +
  geom_col()+
  coord_flip() +
  scale_fill_manual(values = c('Up'= '#36638a','NoSignifi'='#cccccc','Down'='#7bcd7b')) +
  geom_hline(yintercept = c(-1,1),color = 'white',size = 0.5,lty='dashed') +
  xlab('') + 
  ylab('t value of GSVA score, S2 VS S1') + #注意坐标轴旋转了
  guides(fill=F)+ # 不显示图例
  theme_prism(border = T) +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
p
# 添加标签

# 小于-1的数量
low1 <- dat_plot %>% filter(t < -21) %>% nrow()
# 小于0总数量
low0 <- dat_plot %>% filter( t < 0) %>% nrow()
# 小于1总数量
high0 <- dat_plot %>% filter(t < 1) %>% nrow()
# 总的柱子数量
high1 <- nrow(dat_plot)

# 依次从下到上添加标签
p <- p + geom_text(data = dat_plot[1:low1,],aes(x = id,y = 0.1,label = id),
                      hjust = 0,color = 'black') + # 小于-1的为黑色标签
  geom_text(data = dat_plot[(low1 +1):low0,],aes(x = id,y = 0.1,label = id),
            hjust = 0,color = 'grey') + # 灰色标签
  geom_text(data = dat_plot[(low0 + 1):high0,],aes(x = id,y = -0.1,label = id),
            hjust = 1,color = 'grey') + # 灰色标签
  geom_text(data = dat_plot[(high0 +1):high1,],aes(x = id,y = -0.1,label = id),
            hjust = 1,color = 'black') # 大于1的为黑色标签

p

ggsave("gsva_bar.pdf",p,width = 8,height  = 8)




#######################################################################################

## GO富集分析

## 找差异基因

# 1. 定义处理组和对照组的标识（这里定义了一个包含“处理组样本标识的向量”和一个包含“对照组样本标识的向量”。）  
treatments <- c("Human_Organoid_AD","Human_3D_AD", "Human_AD_01","Human_AD_02","Human_AD_03","Human_AD_04","Human_AD_05","Human_AD_06",
                "Human_AD_07","Human_AD_08","Human_AD_09","Human_AD_10")  
controls <- c("Human_Organoid_Normal","Human_3D_Normal", "Human_Normal_01","Human_Normal_02","Human_Normal_03","Human_Normal_04","Human_Normal_05",
              "Human_Normal_06","Human_Normal_07")   

# 2. 创建一个新的标识列，将处理组合并为一个组，对照组合并为一个组  
scRNA_harmony[["treat_vs_control"]] <- 
  ifelse(scRNA_harmony$orig.ident %in% treatments, #判断每个样本的原始标识（orig.ident）是否在处理组的标识中
         "Treated",
         "Control")  

# 3. 将"treat_vs_control"设置为astrocyte_subset新的标识
Idents(scRNA_harmony) <- "treat_vs_control"

# 4. 只保留（提取）标识为"Treated"或"Control"的样本
scRNA_harmony <- subset(scRNA_harmony, idents = c("Treated", "Control"))

# 5. 通过JoinLayers函数合并不同层的数据
scRNA_harmony <- JoinLayers(scRNA_harmony)

# 6. 使用新的标识列"treat_vs_control"进行差异表达分析  
degdf <- FindMarkers(scRNA_harmony, ident.1 = "Treated",   
                                    ident.2 = "Control",   
                                    min.pct = 0.1,   
                                    logfc.threshold = 0.25)  
table(scRNA_harmony$orig.ident)
table(scRNA_harmony$celltype, scRNA_harmony$orig.ident)

# table(scRNA_harmony$orig.ident)
# table(scRNA_harmony$celltype,scRNA_harmony$orig.ident)
# degdf <-  FindMarkers(scRNA_harmony,ident.1 = "sample2",ident.2 = "sample21", 
#                         logfc.threshold = 0.5,group.by = "orig.ident",
#                         ident=1,subset.ident = "Chondrocytes",min.pct = 0.4)

## 筛选p值不显著的基因
degdf1 = filter(degdf, p_val_adj < 0.01)

## options(download.file.method = "curl") 查看当前下载进度条
##BiocManager::install("org.Hs.eg.db")

## GO富集分析
library(org.Hs.eg.db)
library(clusterProfiler)
degs.list = rownames(degdf)
erich.go.BP = enrichGO(gene =degs.list,
                       OrgDb = org.Hs.eg.db,
                       keyType = "SYMBOL",
                       ont = "BP",
                       pvalueCutoff = 0.05,
                       qvalueCutoff = 0.05)

## 气泡图
p1 <- dotplot(erich.go.BP,showCategory = 10 )

ggsave("erich_go_BP.pdf", p1, width = 8, height  = 8)
?dotplot

## barplot
p2 <- barplot(erich.go.BP,showCategory = 8)

ggsave("erich_go_BP_barplot.pdf", p2, width = 8, height  = 8)

erich.go.BP=erich.go.BP@result

setwd("/home/oelab/RunData2/heisenberg/AllDataResult/")
write.table(erich.go.BP,"erich.go.BP.deg.con.hf.txt",sep = "\t",col.names = NA)


######################################################################
library(msigdbr)
msigdbr_species()
msigdbr_collections()
##pathway=read.gmt("c2.cp.kegg.v7.5.1.symbols.gmt")
library(msigdbr)
m_df.h <- msigdbr(species ="Homo sapiens",  category = "H" )

## 提取m_df.h 的第三列和第四列的数据
m_df.h = m_df.h[,c(3,4)] 

?enricher
res = enricher(degs.list, TERM2GENE =m_df.h )

dotplot(res)


### kegg
m_df.kegg <- msigdbr(species ="Homo sapiens",  category = "C2", subcategory = "KEGG" )
m_df.kegg = m_df.kegg[,c(3,4)] 


## enricher 既可以计算GO，也可以计算KEGG，也可以计算HALLMARK
res.kegg = enricher( degs.list,TERM2GENE =m_df.kegg )

dotplot( res.kegg)


############################################################################################################################################
############################################################################################################################################
############################################################################################################################################

######################################################################################################

## 富集分析结果用合旋图展示

colnames(erich.go.BP)
### install.packages("GOplot")
library(GOplot)

## > colnames(erich.go.BP)
## [1] "ID"             "Description"    "GeneRatio"      "BgRatio"        "RichFactor"    
## [6] "FoldEnrichment" "zScore"         "pvalue"         "p.adjust"       "qvalue"        
## [11] "geneID"         "Count"       
go1 = erich.go.BP[1:3,c(1,2,11,9)]
rownames(go1)=go1$ID
go = go1
###install.packages("stringr")
library(stringr)
## 把 / 替换为 ,
go$geneID=str_replace_all(go$geneID,"/",",")
names(go)=c('ID','Term','Genes','adj_pval')
go$Category="BP"



x1=strsplit(go$Genes[1],split=",",fixed=T)
x2=strsplit(go$Genes[2],split=",",fixed=T)
x3=strsplit(go$Genes[3],split=",",fixed=T)
g1=c(x1[[1]],x2[[1]],x3[[1]])


genedata1=degdf[g1,]   
genedata1$ID=rownames(genedata1)
genedata2=data.frame(ID=genedata1$ID,logFC=genedata1$avg_log2FC)
#genedata=data.frame(ID=degs.list,logFC=degdf[degs.list,]$avg_log2FC)


circ <- circle_dat(go,genedata2)
##条形图
GOBar(subset(circ, category == 'BP'))
#气泡图
GOBubble(circ, labels = 3)

GOCircle(circ, nsub = 3)


chord <- chord_dat(circ, genedata2)

chord

GOChord(chord, gene.order = 'logFC')

##基因与GO Term的热图(GOHeat)
GOHeat(chord, nlfc = 1, fill.col = c('red', 'yellow', 'green'))




#################################################################3
############################################################################################################################################
############################################################################################################################################
############################################################################################################################################
############################################################################################################################################
############################################################################################################################################
############################################################################################################################################
############################################################################################################################################
############################################################################################################################################
############################################################################################################################################


library(org.Mm.eg.db)
library(clusterProfiler)
degs.list=rownames(degdf)
erich.go.BP = enrichGO(gene =degs.list,
                       OrgDb = org.Mm.eg.db,
                       keyType = "SYMBOL",
                       ont = "BP",
                       pvalueCutoff = 0.05,
                       qvalueCutoff = 0.05)
dotplot(erich.go.BP,showCategory = 10 )
?dotplot
barplot(erich.go.BP,showCategory = 8)
erich.go.BP=erich.go.BP@result
write.table(erich.go.BP,"4.7.erich.go.BP.deg.con.hf.txt",sep = "\t",col.names = NA)




kegg=read.table("4.7.go.txt",header = T,sep = "\t")
View(kegg)

k = data.frame(kegg)
library(ggplot2)
library(dplyr)
before &lt;- as.numeric(sub("/\\d+$", "", k$GeneRatio))
after &lt;- as.numeric(sub("^\\d+/", "", k$GeneRatio))
k$GeneRatio = before /after
font.size =10

k %&gt;% 
  ## 对进行p值排序
  arrange(p.adjust) %&gt;% 
  ##指定富集的通路数目
  dplyr::slice(1:14) %&gt;% 
  ## 开始ggplot2 作图，其中fct_reorder调整因子level的顺序
  ggplot(aes(GeneRatio,forcats::fct_reorder(Description,Count)))+ 
  ## 画出点图
  geom_point(aes(color=p.adjust, size = Count)) +
  ## 调整颜色，guide_colorbar调整色图的方向
  scale_color_continuous(low="red", high="blue", guide=guide_colorbar(reverse=TRUE))+
  ## 调整泡泡的大小
  scale_size_continuous(range=c(3, 8))+
  ## 如果用ylab("")或出现左侧空白
  labs(y=NULL) +
  ## 如果没有这一句，上方会到顶
  ggtitle("")+
  ## 设定主题
  theme_bw() +
  theme(axis.text.x = element_text(colour = "black",
                                   size = font.size, vjust =1 ),
        axis.text.y = element_text(colour = "black",
                                   size = font.size, hjust =1 ),
        axis.title = element_text(margin=margin(10, 5, 0, 0),
                                  color = "black",size = font.size),
        axis.title.y = element_text(angle=90))




erich.go.CC = enrichGO(gene =degs.list,
                       OrgDb = org.Mm.eg.db,
                       keyType = "SYMBOL",
                       ont = "CC",
                       pvalueCutoff = 0.5,
                       qvalueCutoff = 0.5)
dotplot(erich.go.CC,showCategory = 8)
?dotplot
barplot(erich.go.CC,showCategory = 8)

erich.go.MF = enrichGO(gene =degs.list,
                       OrgDb = org.Mm.eg.db,
                       keyType = "SYMBOL",
                       ont = "MF",
                       pvalueCutoff = 0.5,
                       qvalueCutoff = 0.5)
dotplot(erich.go.MF,showCategory = 8)
?dotplot
barplot(erich.go.MF,showCategory = 8)

keytypes(org.Hs.eg.db)


DEG.entrez_id = mapIds(x = org.Mm.eg.db,
                       keys =  degs.list,
                       keytype = "SYMBOL",
                       column = "ENTREZID")

## install.packages("R.utils")

library(R.utils)

R.utils::setOption("clusterProfiler.download.method","auto")

erich.kegg.res &lt;- enrichKEGG(gene = DEG.entrez_id,
                                organism = "mmu",
                                keyType = "kegg")

dotplot(erich.kegg.res)
kegg.res.df=erich.kegg.res@result
write.table(kegg.res.df,"erich.kegg.res.txt",sep = "\t",col.names = NA)

拟时间分析

伪时间分析：Monocle 能够推断出细胞的发育轨迹，也就是在没有时间点信息的前提下，基于基因表达变化预测细胞随时间的状态演变。这种分析可以帮助理解细胞在发育或分化过程中的变化路径。可以查看 处理组与对照组之间的动态变化轨迹（即：control组到处理组的变化），也可以查看 细胞的发育过程（即：Naive B细胞到 Memory B细胞的变化），也就是 细胞从一种状态到另外一种状态变化的时候，细胞经历了什么过程

Monocle2

Monocle2 是利用机器学习技术（反向图嵌入）从单细胞数据中学习显式的主图来对细胞进行排序（例如从Naive B到 memory B的排序），主要还是用Monocle2进行轨迹的构建

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236

## Monocle2代码

##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())

setwd("/home/oelab/RunData2/heisenberg/")
scRNA_harmony <- readRDS("/home/oelab/RunData2/heisenberg/scRNA_harmony.rds")

library(dplyr)
library(Seurat)
library(tidyverse)
library(patchwork)

table(scRNA_harmony@meta.data$celltype)
table(scRNA_harmony@meta.data$seurat_clusters)

scRNA_harmony <- JoinLayers(scRNA_harmony)

## 单独提取 xx 细胞进行轨迹的分析
## 以Micro.为例
sc.t = scRNA_harmony[,rownames(subset(scRNA_harmony@meta.data,celltype == "Micro."))]  

options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")

# BiocManager::install("monocle")
library(monocle)
Idents(scRNA_harmony)="celltype"

## 这里的拟时序分析可以做对照组到处理组的，也可以做两个细胞之间的变化
scRNA.Micro = subset(scRNA_harmony,ident=c("Micro."))

table(scRNA.Micro$orig.ident)

# 定义处理组和对照组的标识（这里定义了一个包含“处理组样本标识的向量”和一个包含“对照组样本标识的向量”。）  
# treatments <- c("Human_Organoid_AD","Human_3D_AD", "Human_AD_01","Human_AD_02","Human_AD_03","Human_AD_04","Human_AD_05","Human_AD_06",
#                 "Human_AD_07","Human_AD_08","Human_AD_09","Human_AD_10")  
# controls <- c("Human_Organoid_Normal","Human_3D_Normal", "Human_Normal_01","Human_Normal_02","Human_Normal_03","Human_Normal_04","Human_Normal_05",
#               "Human_Normal_06","Human_Normal_07")   
# 
# 
# # 创建一个新的标识列，将处理组合并为一个组，对照组合并为一个组  
# scRNA.Micro[["treat_vs_control"]] <- 
#   ifelse(scRNA.Micro$orig.ident %in% treatments, #判断每个样本的原始标识（orig.ident）是否在处理组的标识中
#          "Treated",
#          "Control")  
# 
# # 设为新的标识
# Idents(scRNA.Micro) <- "treat_vs_control"
# 
# ##
# ## Idents(scRNA.Micro)="orig.ident"
# sc.Micro = subset(scRNA.Micro, ident=c("Treated", "Control"))
# table(sc.Micro$orig.ident)

######
## past0是连接符
## 提取对应的细胞的流程
# scRNA_harmony$huage=paste0(scRNA_harmony$orig.ident,"_",scRNA_harmony$celltype)
# table(scRNA_harmony$huage)
# Idents(scRNA_harmony)="huage"
# sc.123=subset(scRNA.Osteoclastic,ident="sample2_Chondrocytes")

## 查看可用的层次
scRNA.Micro[["RNA"]]@layers

data <- GetAssayData(scRNA.Micro,slot = "count")

data[1:20,1:20]

?new
## new就是转换成S4对象
pd <- new('AnnotatedDataFrame', data = scRNA.Micro@meta.data)
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)
monocle_cds <- newCellDataSet(data,
                                 phenoData = pd,
                                 featureData = fd,
                                 lowerDetectionLimit = 0.5,
                                 expressionFamily = negbinomial.size())



monocle_cds <- estimateSizeFactors(monocle_cds)
monocle_cds <- estimateDispersions(monocle_cds)

monocle_cds <- detectGenes(monocle_cds, min_expr = 0.1)
print(head(fData(monocle_cds)))

HSMM = monocle_cds
## Monocle2高变基因的筛选
disp_table <- dispersionTable(HSMM)
disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit)$gene_id
HSMM <- setOrderingFilter(HSMM, disp.genes)
plot_ordering_genes(HSMM)

## 降维
## max_components = 2 指两个维度
HSMM <- reduceDimension(HSMM, max_components = 2,
                           method = 'DDRTree')

## 排序，指明发育起点和发育终点
HSMM <- orderCells(HSMM)

## 图示说明
plot_cell_trajectory(HSMM, color_by = "seurat_clusters")

plot_cell_trajectory(HSMM, color_by = "celltype")

plot_cell_trajectory(HSMM, color_by = "orig.ident")

plot_cell_trajectory(HSMM, color_by = "State")

plot_cell_trajectory(HSMM, color_by = "Pseudotime")

# root_state = 5：这一步是指定拟时分析的“根状态”，也就是定义拟时过程的起点。
# 通过设置 root_state = 5，你指定第 5 个状态为细胞分化或发育过程的起始点。
# 在没有指定 root_state 时，Monocle 会自动选择一个根状态
HSMM <- orderCells(HSMM, root_state = 5)

plot_cell_trajectory(HSMM, color_by = "Pseudotime")



plot_cell_trajectory(HSMM, color_by = "celltype")

## facet_wrap 用来做分支图
## 展示不同clusters的状态
plot_cell_trajectory(HSMM, color_by = "State") +
  facet_wrap(~State, nrow = 3)

plot_cell_trajectory(HSMM, color_by = "State") +
  facet_wrap(~orig.ident )

## 展示基因的变化的散点图
## "GAPDH", "RORA"为示例基因
blast_genes <- row.names(subset(fData(HSMM),
                                   gene_short_name %in% c("GAPDH", "RORA")))

plot_genes_jitter(HSMM[blast_genes,],
                  grouping = "State",
                  min_expr = 0.1)

## 基因随着发育时间的变化图 
## 筛选基因
HSMM_expressed_genes <-  row.names(subset(fData(HSMM),
                                             num_cells_expressed >= 10))

HSMM_filtered <- HSMM[HSMM_expressed_genes,]

my_genes <- row.names(subset(fData(HSMM_filtered),
                                gene_short_name %in% c("YWHAB", "GAPDH", "TNNC1")))
cds_subset <- HSMM_filtered[my_genes,]
plot_genes_in_pseudotime(cds_subset, color_by = "State")

plot_genes_in_pseudotime(cds_subset, color_by =  "celltype")


## 想看下列三个基因的小提琴图、pseudotime图等
genes <- c("TNNT2", "TNNC1", "CDK1")
p1 <- plot_genes_jitter(HSMM[genes,], grouping = "State", color_by = "State")
p2 <- plot_genes_violin(HSMM[genes,], grouping = "State", color_by = "State")
p3 <- plot_genes_in_pseudotime(HSMM[genes,], color_by = "State")
plotc <- p1|p2|p3
plotc 



to_be_tested <- row.names(subset(fData(HSMM),
                                    gene_short_name %in% c("MYH3", "MEF2C", "CCNB2", "TNNT1")))
cds_subset <- HSMM[to_be_tested,]

diff_test_res <- differentialGeneTest(cds_subset,
                                         fullModelFormulaStr = "~sm.ns(Pseudotime)")

diff_test_res[,c("gene_short_name", "pval", "qval")]
plot_genes_in_pseudotime(cds_subset, color_by ="State")



## 筛选下列的基因的拟时序分析
marker_genes <- row.names(subset(fData(HSMM),
                                    gene_short_name %in% c("MEF2C", "MEF2D", "MYF5",
                                                           "ANPEP", "PDGFRA","MYOG",
                                                           "TPM1",  "TPM2",  "MYH2",
                                                           "MYH3",  "NCAM1", "TNNT1",
                                                           "TNNT2", "TNNC1", "CDK1",
                                                           "CDK2",  "CCNB1", "CCNB2",
                                                           "CCND1", "CCNA1", "ID1")))

diff_test_res <- differentialGeneTest(HSMM[marker_genes,],
                                         fullModelFormulaStr = "~sm.ns(Pseudotime)")

sig_gene_names <- row.names(subset(diff_test_res, qval < 1))

plot_pseudotime_heatmap(HSMM[sig_gene_names,],
                        num_clusters = 4,
                        cores = 1,
                        show_rownames = T)


plot_cell_trajectory(HSMM, color_by = "State")

################################################################################

## 无法运行！！！
## 查看哪些基因向左升高，哪些基因向右升高
BEAM_res <- BEAM(HSMM, branch_point = 1, progenitor_method = "duplicate", cores = 7)

# methods <- c("duplicate", "expression", "cluster")
# results <- lapply(methods, function(method) {
#   BEAM(HSMM, branch_point = 1, progenitor_method = method)
# })

BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
plot_genes_branched_heatmap(HSMM[row.names(subset(BEAM_res,
                                                  qval < 1e-4)),],
                            branch_point = 1,
                            num_clusters = 4,
                            cores = 1,
                            use_gene_short_name = T,
                            show_rownames = T)


genes <- row.names(subset(fData(HSMM),
                             gene_short_name %in% c( "MEF2C", "CCNB2", "TNNT1")))

plot_genes_branched_pseudotime(HSMM[genes,],
                               branch_point = 1,
                               color_by = "State",
                               ncol = 1)

Monocle3

与Monocle2相比，Monocle3最大的变化就是可以去除批次效应，并且Monocle2有细胞数量的限制，Monocle3具有环路或收敛点的轨迹的方法

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286

## Monocle3的代码

## BiocManager::install(c('BiocGenerics', 'DelayedArray', 'DelayedMatrixStats',
## 'limma', 'S4Vectors', 'SingleCellExperiment',
## 'SummarizedExperiment', 'batchelor', 'Matrix.utils'))



## install.packages("devtools")
## devtools::install_github('cole-trapnell-lab/leidenbase')
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")

options(download.file.method = "curl")

## devtools::install_github('cole-trapnell-lab/monocle3')


################################################################################


##系统报错改为英文
Sys.setenv(LANGUAGE = "en")
##禁止转化为因子
options(stringsAsFactors = FALSE)
##清空环境
rm(list=ls())

library(dplyr)
library(Seurat)
library(tidyverse)
library(patchwork)
library(monocle3)

setwd("/home/oelab/RunData2/heisenberg/")

scRNA_harmony <- readRDS("/home/oelab/RunData2/heisenberg/scRNA_harmony.rds")
scRNA_harmony <- JoinLayers(scRNA_harmony)


setwd("/home/oelab/RunData2/heisenberg/AllDataResult/")

data <- GetAssayData(scRNA_harmony, slot = 'counts')
cell_metadata <- scRNA_harmony@meta.data
gene_annotation <- data.frame(gene_short_name = rownames(data))
rownames(gene_annotation) <- rownames(data)
cds <- new_cell_data_set(data,
                            cell_metadata = cell_metadata,
                            gene_metadata = gene_annotation)

## 归一化处理
## 这个时间比较久
cds <- preprocess_cds(cds, num_dim = 100)

## 展示前100个pc
plot_pc_variance_explained(cds)

? reduce_dimension

## 降维 umap
cds <- reduce_dimension(cds)

## 降维画图
plot_cells(cds)
DimPlot(scRNA_harmony)
plot_cells(cds, color_cells_by="seurat_clusters")

## 查看某个基因的表达量
plot_cells(cds, genes=c("GAPDH"))

?reduce_dimension
## 进行tsne降维
cds <- reduce_dimension(cds,umap.fast_sgd=T, cores = 3, reduction_method="tSNE")

plot_cells(cds, reduction_method="tSNE")

plot_cells(cds, reduction_method="tSNE", color_cells_by="celltype")

plot_cells(cds, color_cells_by="orig.ident", label_cell_groups=FALSE, reduction_method="tSNE")

## alignment_group指明消除批次效应的来源
cds = align_cds(cds, num_dim = 100, alignment_group ="orig.ident")
cds = reduce_dimension(cds, reduction_method="tSNE")
plot_cells(cds, color_cells_by="orig.ident", label_cell_groups=FALSE,,reduction_method="tSNE")

cds = reduce_dimension(cds)

## 进行clusters聚类
cds = cluster_cells(cds, resolution=1e-5)
plot_cells(cds)

table(cds@clusters$UMAP$clusters)
cds@clusters$UMAP$partitions
plot_cells(cds, color_cells_by="partition", group_cells_by="partition")

plot_cells(cds, color_cells_by="celltype")

marker_test_res <- top_markers(cds, group_cells_by="partition", 
                                  reference_cells=1000, cores=8)


top_specific_markers <-  marker_test_res %>%
  filter(fraction_expressing >= 0.10) %>%
  group_by(cell_group) %>%
  top_n(1, pseudo_R2)

top_specific_marker_ids <- unique(top_specific_markers %>% pull(gene_id))

plot_genes_by_group(cds,
                    top_specific_marker_ids,
                    group_cells_by="celltype",
                    ordering_type="maximal_on_diag",
                    max.size=3)





top_specific_markers = marker_test_res %>%
  filter(fraction_expressing >= 0.10) %>%
  group_by(cell_group) %>%
  top_n(3, pseudo_R2)

top_specific_marker_ids = unique(top_specific_markers %>% pull(gene_id))

plot_genes_by_group(cds,
                    top_specific_marker_ids,
                    group_cells_by="celltype",
                    ordering_type="cluster_row_col",
                    max.size=3)

# 先将partitions的分组由因子型转为字符型
colData(cds)$assigned_cell_type <- as.character(partitions(cds))

cds_subset <- choose_cells(cds)


cds <- learn_graph(cds)
plot_cells(cds,
           color_cells_by = "celltype",
           label_groups_by_cluster=FALSE,
           label_leaves=FALSE,
           label_branch_points=FALSE,
           group_label_size=4,cell_size=0.5)


cds = order_cells(cds)

plot_cells(cds,
           color_cells_by = "pseudotime",
           label_cell_groups=FALSE,
           label_leaves=TRUE,
           label_branch_points=F,
           graph_label_size=1.5,
           group_label_size=4,cell_size=0.5)


table(scRNA_harmony@meta.data$celltype)
# a helper function to identify the root principal points:
get_earliest_principal_node <- function(cds, time_bin="T_cells"){
  cell_ids <- which(colData(cds)[, "celltype"] == time_bin)
  
  closest_vertex <- cds@principal_graph_aux[["UMAP"]]$pr_graph_cell_proj_closest_vertex
  closest_vertex <- as.matrix(closest_vertex[colnames(cds), ])
  root_pr_nodes <-
    igraph::V(principal_graph(cds)[["UMAP"]])$name[as.numeric(names(which.max(table(closest_vertex[cell_ids,]))))]
  
  root_pr_nodes
}
cds = order_cells(cds, root_pr_nodes=get_earliest_principal_node(cds))

plot_cells(cds,
           color_cells_by = "pseudotime",
           label_cell_groups=FALSE,
           label_leaves=FALSE,
           label_branch_points=FALSE,
           graph_label_size=1.5,
           group_label_size=4,cell_size=1.5)


cds_sub <- choose_graph_segments(cds)




# Working with 3D trajectories
## 制作三维的Monocle3图

cds_3d = reduce_dimension(cds, max_components = 3)
cds_3d = cluster_cells(cds_3d)
cds_3d = learn_graph(cds_3d)
cds_3d = order_cells(cds_3d, root_pr_nodes=get_earliest_principal_node(cds))

cds_3d = order_cells(cds_3d)

cds_3d_plot_obj = plot_cells_3d(cds_3d, color_cells_by="celltype")
cds_3d_plot_obj


ciliated_genes = top_specific_markers$gene_id[5:10]
cds_subset = cds[rowData(cds)$gene_short_name %in% ciliated_genes,]

gene_fits = fit_models(cds_subset, model_formula_str = "~seurat_clusters")
fit_coefs
fit_coefs = coefficient_table(gene_fits)
# 挑出时间相关的组分
emb_time_terms = fit_coefs %>% filter(term == "seurat_clusters1")
emb_time_terms

emb_time_terms = fit_coefs %>% filter(term == "seurat_clusters1")

emb_time_terms %>% filter (q_value < 0.05) %>%
  select(gene_short_name, term, q_value, estimate)

plot_genes_violin(cds_subset[,], group_cells_by="celltype", ncol=2) +
  theme(axis.text.x=element_text(angle=45, hjust=1))


neurons_cds <- cds[,grepl("Macrophage", colData(cds)$celltype, ignore.case=TRUE)]
plot_cells(neurons_cds, color_cells_by="partition")

pr_graph_test_res <- graph_test(neurons_cds, neighbor_graph="knn", cores=3)
pr_deg_ids <- row.names(subset(pr_graph_test_res, q_value < 0.05))

gene_module_df = find_gene_modules(neurons_cds[pr_deg_ids,], resolution=1e-2)

cell_group_df = tibble::tibble(cell=row.names(colData(neurons_cds)), cell_group=neurons_cds@colData@listData[["seurat_clusters"]])
agg_mat = aggregate_gene_expression(neurons_cds, gene_module_df, cell_group_df)
row.names(agg_mat) = stringr::str_c("Module ", row.names(agg_mat))
colnames(agg_mat) = stringr::str_c("Partition ", colnames(agg_mat))

pheatmap::pheatmap(agg_mat, cluster_rows=TRUE, cluster_cols=TRUE,
                   scale="column", clustering_method="ward.D2",
                   fontsize=6)


plot_cells(neurons_cds,
           genes=gene_module_df %>% filter(module %in% c(16,38,33,42)),
           group_cells_by="cluster",
           color_cells_by="cluster",
           show_trajectory_graph=FALSE)


plot_cells(cds,
           color_cells_by = "celltype",
           label_groups_by_cluster=FALSE,
           label_leaves=FALSE,
           label_branch_points=FALSE)


##找到影响发育轨迹的基因
trace('calculateLW', edit = T, where = asNamespace("monocle3"))
## change Matrix::rBind to rbind on line 93
# 使用neighbor_graph="principal_graph"来检验轨迹相邻的细胞的表达是否相关
ciliated_cds_pr_test_res <- graph_test(cds, neighbor_graph="principal_graph", cores=4)



pr_deg_ids <- row.names(subset(ciliated_cds_pr_test_res, q_value < 0.05))

AFD_genes <- pr_deg_ids[20:22]
AFD_lineage_cds <- cds[rowData(cds)$gene_short_name %in% AFD_genes,
]
plot_genes_in_pseudotime(AFD_lineage_cds,
                         color_cells_by="celltype",
                         min_expr=0.5)


cds_subset <- choose_cells(cds)

subset_pr_test_res <- graph_test(cds_subset, neighbor_graph="principal_graph", cores=4)
pr_deg_ids <- row.names(subset(subset_pr_test_res, q_value < 0.05))

gene_module_df <- find_gene_modules(cds_subset[pr_deg_ids,], resolution=0.001)

agg_mat <- aggregate_gene_expression(cds_subset, gene_module_df)
module_dendro <- hclust(dist(agg_mat))
gene_module_df$module <- factor(gene_module_df$module, 
                                   levels = row.names(agg_mat)[module_dendro$order])

plot_cells(cds_subset,
           genes=gene_module_df,
           label_cell_groups=FALSE,
           show_trajectory_graph=FALSE)


细胞互作

CellChat