添加 10xscRNA-seq.R

Signed-off-by: 生信分析 <bioinfo@baihub.cn>
2024-12-11 02:53:16 +08:00 · 2024-12-11 02:53:16 +08:00 · b88b661378
commit b88b661378
parent 2f3e8397e9
1 changed files with 114 additions and 0 deletions
--- a/10xscRNA-seq.R
+++ b/10xscRNA-seq.R
@ -0,0 +1,114 @@
+## 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"))
+