ESCC_LDRT_CI/scRNA_quality.R at main · GuoBioinfoLab/ESCC_LDRT_CI · 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
rm(list=ls())
library(Seurat)
library(ggplot2)
library(tidyverse)
setwd("2-harmony/")
sce.all.int <-  readRDS("sce.all_int.rds")

feats <- c("nFeature_RNA", "nCount_RNA", "percent_mito")
# p1_filtered=VlnPlot(sce.all.int, group.by = "orig.ident", features = feats, pt.size = 0, ncol = 2) +
#   NoLegend()
# w=length(unique(sce.all.int$orig.ident))/3+8;w
# ggsave(filename="Vlnplot1_filtered.pdf",plot=p1_filtered,width = w,height = 5)
meta1 <- sce.all.int@meta.data %>% select(sample,nFeature_RNA,nCount_RNA,percent_mito)
# 宽格式转长格式
meta2 <- meta1 %>%
  pivot_longer(
    cols = c(nFeature_RNA, nCount_RNA, percent_mito),
    names_to = "feature",
    values_to = "value"
  )
#meta2$sample <-  factor(meta2$sample, levels = c(order1))
# 绘制小提琴图
meta <- read.delim("../../../meta.txt")
meta <- meta[meta$class1 == "scRNA",]

meta2$group <- meta$group1[match(meta2$sample,meta$id)]
order1 <- c("Pre_NR","Post_NR","Pre_R","Post_R")
meta2$group <- factor(meta2$group, levels = order1)
color1 <- c("#4682B4","#87CEFA","#BC8F8F","#D2B48C")
p1 <- ggplot(meta2, aes(x = sample, y = value, fill = group)) +
  geom_violin(scale = "width", trim = TRUE) +
  facet_wrap(~ feature, ncol = 1, scales = "free") +
  scale_fill_manual(values = color1)+
  theme_classic() +
  theme(
    strip.background = element_blank(),  # 移除背景框
    strip.text = element_text(size = 14),  # 增大字体大小
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid = element_blank(),      # 删除所有网格线
    panel.background = element_blank() # 删除背景
  ) +
  labs(x = "", y = "")
p1
w=length(unique(sce.all.int$orig.ident))/3
ggsave(p1, filename="Vlnplot1_filtered.pdf",width = w,height = 6)
ggsave(p1, filename="Vlnplot1_filtered.png",width = w,height = 6)


p1 <- DimPlot(sce.all.int, reduction = "umap", group.by = "orig.ident", pt.size = 1) +
  labs(title = "")
p1
ggsave(p1, file = "umap_sample.pdf", height = 5, width = 8)
ggsave(p1, file = "umap_sample.png", height = 5, width = 8)

meta <- read.delim("../../../meta.txt")
meta <- meta[meta$class1 == "scRNA",]
sce.all.int$group <- meta$group1[match(sce.all.int$orig.ident,meta$id)]

sce.all.int$group <- factor(sce.all.int$group, levels = rev(c("Pre_NR","Post_NR","Pre_R","Post_R")))
p1 <- DimPlot(sce.all.int, reduction = "umap", group.by = "group", pt.size = 1) +
  labs(title = "") + theme_void() +
  scale_color_manual(values = c(
    "Pre_NR" = "#4682B4",
    "Post_NR" = "#87CEFA",
    "Pre_R" = "#BC8F8F",
    "Post_R" = "#D2B48C"
  )) +
  guides(color = guide_legend(reverse = TRUE))

p1

ggsave(p1, file = "umap_group.pdf", height = 5, width = 6)
ggsave(p1, file = "umap_group.png", height = 5, width = 6)


twothosuand <- c("S27_C_Pre", "S15_C_Post") #显示特殊样本
sce.all.int$color <- sce.all.int$orig.ident
sce.all.int$color[!(sce.all.int$color %in% twothosuand)] <- "others"
p1 <- DimPlot(sce.all.int, reduction = "umap", group.by = "color", pt.size = 0.5)
#ggsave(p1, file = "umapsmall3.pdf", height = 5, width = 7)

avg <- AverageExpression(object = sce.all.int, group.by = "orig.ident")
avg <- avg$RNA

library(factoextra)
avg <- t(avg)
avg <- as.data.frame(avg)
rownames(avg) #samples are row
avg <- avg[,colSums(avg) > 0]
avg.scaled <- scale(avg) # 对数据进行标准化处理，这是PCA之前的常见步骤

# 执行PCA
pca_res <- prcomp(avg.scaled, center = TRUE, scale. = TRUE)
site <- as.data.frame(pca_res$x[,1:2])
site$id <- rownames(site)

library(ggrepel)
p <- ggplot(data = site, aes(x = PC1, y = PC2)) +
  geom_point() +  #添加样本点
  theme_classic() +
  labs(x = paste0("PC1 (", round(100*summary(pca_res)$importance[2,1], 1), "%)"),
       y = paste0("PC2 (", round(100*summary(pca_res)$importance[2,2], 1), "%)"),) +
  geom_text_repel(aes(label = id),  # 使用 geom_text_repel 避免标签重叠
                  box.padding = 0.35,  # 标签周围的填充
                  point.padding = 0.5, # 标签与点的间距
                  size = 3,            # 标签字体大小
                  color = "blue",      # 标签颜色
                  fontface = "bold") # 标签字体样式


p
ggsave(p, file = "pca.pdf", height = 5, width = 6)

# plot(pca_res$x[,1], pca_res$x[,2],
#      xlab = paste0("PC1 (", round(100*summary(pca_res)$importance[2,1], 1), "%)"),
#      ylab = paste0("PC2 (", round(100*summary(pca_res)$importance[2,2], 1), "%)"),
#      main = "PCA Score Plot")