## Written By LiShang
## Notice: 仅适用于10x平台

#加载依赖包
library(Seurat)
library(sctransform)
library(dplyr)
library(ggplot2)

########################读取样本并生成Seurat Object##########################

#获取所有样本的存放目录
setwd("~/scRNA-seq")
sample_dirs <- list.dirs(full.names = FALSE)[-1]

#同时读取所有样本, 并转化成Seurat Object
sample_list <- mapply(function(dir) {
  Read10X(data.dir = dir) %>%
    CreateSeuratObject(project = dir, min.cells = 3, min.features = 200)
}, dir = sample_dirs, SIMPLIFY = FALSE)
#合并所有样本
merged <- merge(sample_list[[1]], sample_list[-1], add.cell.ids = sample_dirs)

#提取样本的分组信息, 本模板分为cancer和normal两个组
merged[["sample.grp"]] <- merged@meta.data$orig.ident %>%
  strsplit(split = "_") %>%
  sapply(function(x) x[2]) %>%
  as.factor()

#################################质量控制####################################

#计算线粒体转录本占比
merged[["percent.mt"]] <- PercentageFeatureSet(merged, pattern = "^MT-")
#细胞过滤
merged <- merged %>%
  subset(subset = nFeature_RNA > 200 & nFeature_RNA < 2500) %>%
  subset(subset = percent.mt < 5)

###############################标准分析流程##################################

#SCTransform标准化、线性降维
merged <- merged %>%
  SCTransform(vars.to.regress = "percent.mt", verbose = FALSE) %>%
  RunPCA(verbose = FALSE)

#多样本数据整合, harmony去批次
merged <- IntegrateLayers(
  object = merged,
  method = HarmonyIntegration,
  normalization.method = "SCT",
  orig.reduction = "pca",
  new.reduction = "harmony",
  verbose = FALSE
)
merged[["RNA"]] <- JoinLayers(merged[["RNA"]])

#聚类、非线性降维
merged <- merged %>%
  FindNeighbors(dims = 1:30, reduction = "harmony", verbose = FALSE) %>%
  FindClusters(resolution = 2, verbose = FALSE) %>%
  RunUMAP(dims = 1:30, reduction = "harmony", verbose = FALSE)

###############################细胞类型注释##################################

DimPlot(merged, reduction = "umap", label = TRUE)
DotPlot(merged, c("EPCAM", "CDH1", "GATA3", "KRT18", #Epithelial
                  "ESR1", "PGR", "ERBB2", #Breast Cancer
                  "PDGFRA", "COL1A1", "FAP", "DCN", "VIM", #Fibroblast
                  "CD34", "VWF", "CDH5", "ENG", "PECAM1", #Endothelial
                  "PTPRC", #CD45: Immune
                  "CD2", "CD3D", "CD3E", "CD3G",  #T cell
                  "MS4A1", "CD79A", "CD79B", #B cell
                  "KLRF1", "KLRD1", "NKG7", "GNLY" #NK cell
                  ), cols = c("yellow","blue")
) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

new_cluster_ids <- c(
  "T cell",
  "malignant",
  "T cell",
  "T cell",
  "malignant",
  "malignant",
  "malignant",
  "malignant",
  "fibroblast",
  "malignant",
  "fibroblast",
  "malignant",
  "malignant",
  "malignant",
  "malignant",
  "endothelial",
  "fibroblast",
  "T cell",
  "fibroblast",
  "malignant",
  "malignant",
  "malignant",
  "malignant",
  "B cell"
)
names(new_cluster_ids) <- levels(merged)
merged <- RenameIdents(merged, new_cluster_ids)

###############################数据分析作图##################################

merged <- PrepSCTFindMarkers(merged)
markers <- FindMarkers(merged, ident.1 = "malignant")

FeaturePlot(merged, features = c("CREB1","CREB3"))
DotPlot(merged, features = c("CREB1","CREB3"))
VlnPlot(merged, features = c("CREB1"))