# install.packages('R.utils')
rm(list=ls())
options(stringsAsFactors = F) 
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(data.table)
library(stringr)
#install.packages("textshape")
library(textshape)
#setwd("../../")
# Получить список всех файлов rds
file_list <- list.files("./GSE234933_raw/rds/", pattern = ".rds")

# Создайте пустой список для хранилищаSeuratобъектов.
seurat_list <- list()

# Цикл для чтения данных каждого файла rds и создания Seuratобъекта.
for (file in file_list) {
  # Путь к файлу сращивания
  data.path <- paste0("./GSE234933_raw/rds/", file)
  # Чтение данных файлов RDS
  seurat_data <- readRDS(data.path)
  # Создайте Seuratобъект и укажите имя проекта в качестве имени файла (удалите суффикс)
  sample_name <- tools::file_path_sans_ext(basename(file))
  seurat_obj <- CreateSeuratObject(counts = seurat_data,
                                   project = sample_name,
                                   min.features = 200,
                                   min.cells = 3)
  # Добавить Seuratобъект в список
  seurat_list <- append(seurat_list, seurat_obj)
}

# Извлеките часть перед подчеркиванием
sample_names <- sub("_.*", "", file_list)
# Объединить объекты Seurat, объединить все объекты Seurat в один объект
seurat_combined <- merge(seurat_list[[1]],
                         y = seurat_list[-1])
# Распечатать объединенный объект Seurat
print(seurat_combined)
seurat_combined@assays$RNA@counts[1:10, 1:2]
##хранилищеданные Это займет больше времени####
saveRDS(seurat_combined, file = "seurat2.rds")

###Прочитайте информацию о субпопуляции ячеек, заданную набором данных.
cell<-fread("./GSE234933_raw/GSE234933_MGH_HNSCC_cell_annotation.txt.gz")
colnames(cell)
meta2 <- cell %>% column_to_rownames("sample_barcode")
###Два метода могут добавлять информацию о подгруппе ячеек, столбец подкластеризации, указанный в статье.
# sce_obj <- CreateSeuratObject(counts = seurat_combined@assays$RNA$counts,
#                               meta.data =meta2,
#                               project = sample_name,
#                               min.features = 200,
#                               min.cells = 3)
sce3 <- AddMetaData(object = seurat_combined, metadata = meta2)
head(sce3@meta.data)

###########
###Группировка и гармония уменьшения размерности
dir.create("2-harmony2")
getwd()
setwd("2-harmony2")

# sce.all=readRDS("../1-QC/sce.all_qc.rds")
sce=sce3 
sce
sce <- NormalizeData(sce, 
                     normalization.method = "LogNormalize",
                     scale.factor = 1e4) 
sce <- FindVariableFeatures(sce)
sce <- ScaleData(sce)
sce <- RunPCA(sce, features = VariableFeatures(object = sce))
#install.packages("harmony")
library(harmony)
seuratObj <- RunHarmony(sce, "orig.ident")
names(seuratObj@reductions)
seuratObj <- RunUMAP(seuratObj,  dims = 1:15, 
                     reduction = "harmony")
DimPlot(seuratObj,reduction = "umap",label=T ) 

sce=seuratObj
sce <- FindNeighbors(sce, reduction = "harmony",
                     dims = 1:15) 
sce.all=sce

sce.all=FindClusters(sce.all, #graph.name = "CCA_snn",
                     resolution = 0.05, algorithm = 1)

sel.clust = "RNA_snn_res.0.05"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident) 
saveRDS(sce.all, "sce.all_int0.05.rds")

####################
#####Проверьте статус группировки
DimPlot(sce.all, reduction = "umap", group.by = "seurat_clusters",label = T)
colnames(sce.all@meta.data)
DimPlot(sce.all, reduction = "umap", group.by = "subclustering",label = T)

##### Подгруппа с удаленной частью NA#######
table(sce.all$subclustering)
###подгруппировать<NA>Назначение0
sce.all$subclustering = sce.all$subclustering %>%
  replace(x = ., list =is.na(.), values =0)
sce.all$celltype<-ifelse(sce.all$subclustering=="Endothelial cells","Endothelial cells",
                              ifelse(sce.all$subclustering=="Fibroblasts","Fibroblasts",
                                     ifelse(sce.all$subclustering=="Lymphocytes","Lymphocytes",
                                            ifelse(sce.all$subclustering=="Myeloid cells","Myeloid cells",
                                                   ifelse(sce.all$subclustering=="Tumor cells","Tumor cells","NA")))))
table(sce.all$celltype)
sce <- sce.all[, !(sce.all$celltype %in% c('NA'))]
DimPlot(sce, reduction = "umap", group.by = "celltype",label=T) 
ggsave('umap_by_RNA_snn_res.0.05_paper.pdf',width =7,height = 6)