scATAC联合scRNA之signac分析(二): 构建scRNA数据及注释-再回顾一下scRNA分析流程吧

setwd('./Documents/KS_TS/signac_scATAC/scRNA_for_inter/')
###################################### load packages#####################################
library(dplyr)library(Seurat)library(patchwork)library(stringr)library(ggplot2)

###################################### process data######################################批量读入data_dirs <- list.dirs(path='.', full.names=TRUE, recursive=FALSE)samples <- c("AA","HC","SD")objs <- list()for(ix in seq_along(data_dirs)){  sample <- samples[ix]  path <- data_dirs[ix]  data <- Read10X(data.dir=path)  obj <- CreateSeuratObject(counts=data, project=sample, min.cells=3, min.features=200) # orig.ident not getting assigned correctly?  obj$orig.ident <- sample  objs[[sample]] <- obj}
#线粒体基因比例for(i in seq_along(objs)){  objs[[i]][["percent.mt"]] <- PercentageFeatureSet(objs[[i]], pattern = "^MT-")}
#filter cellsfor(i in seq_along(objs)){  objs[[i]] <- subset(objs[[i]],                 subset = (                  nFeature_RNA > 200 &                     nFeature_RNA < Inf &                     nCount_RNA > 1000 &                     nCount_RNA < Inf &                     percent.mt < 20                ))}
#run cluster for each sampleseuratPreProcess <- function(obj, selectionMethod="vst", nFeatures=2000, dims=1:10, res=0.5, seed=1){  # perform standard Seurat preprocessing   # (Normalization, Scaling, Dimensionality Reduction, Clustering)  obj <- NormalizeData(obj)  obj <- ScaleData(obj, verbose=FALSE)  obj <- FindVariableFeatures(obj, selection.method=selectionMethod, nfeatures=nFeatures)  obj <- RunPCA(obj)  obj <- RunUMAP(obj, dims=dims)  obj <- FindNeighbors(obj, dims=dims)  obj <- FindClusters(obj, resolution=res, random.seed=1)  return(obj)}

objs <- lapply(seq_along(objs), function(i){    obj <- objs[[i]]  seuratPreProcess(obj)  })

names(objs) <- samples

###################################### Doublet cells detect#####################################source('./ks_detectDoublet')objs[[1]] <- ks_detectDoublet(objs[[1]],dims = 1:15,estDubRate=0.065,                                       ncores = 2, SCTransform=F, Homotypic=F, annotation="seurat_clusters")
objs[[2]] <- ks_detectDoublet(objs[[2]],dims = 1:15,estDubRate=0.025,                                       ncores = 2, SCTransform=F, Homotypic=F, annotation="seurat_clusters")
objs[[3]] <- ks_detectDoublet(objs[[3]],dims = 1:15,estDubRate=0.060,                                       ncores = 2, SCTransform=F, Homotypic=F, annotation="seurat_clusters")

###################################### re-cluster#####################################
# Merge individual Seurat objectsobj <- merge(x=objs[[1]], y=objs[2:length(objs)], add.cell.ids=names(objs))
# Store sample information in metadataobj$Sample <- obj$orig.ident
# Store disease status in metadata# obj$diseaseStatus <- NA# obj$diseaseStatus <- ifelse(grepl("C_SD", obj$Sample), "C_SD", obj$diseaseStatus)# obj$diseaseStatus <- ifelse(grepl("C_PB", obj$Sample), "C_PB", obj$diseaseStatus)# obj$diseaseStatus <- ifelse(grepl("AA", obj$Sample), "AA", obj$diseaseStatus)
# Remove doubletsobj <- subset(obj, subset = (DF.classify == "Singlet"))
#clusterobj <- NormalizeData(obj, normalization.method = "LogNormalize") %>% FindVariableFeatures(., selection.method = "vst", nfeatures = 3000)obj <- ScaleData(obj, vars.to.regress = c("percent.mt"), verbose = T)obj <- RunPCA(obj, npcs = 50, verbose = FALSE)
#run harmonylibrary(harmony)obj <- RunHarmony(obj, "orig.ident")obj <- FindNeighbors(obj, reduction = "harmony",dims = 1:25) obj  <- FindClusters(object = obj , resolution = seq(from = 0.2, to = 1.0,  by = 0.2))obj <- RunUMAP(obj,  dims = 1:25, reduction = "harmony")
# library(clustree)# clustree(obj)

Idents(obj) <- "RNA_snn_res.0.4"obj$seurat_clusters <- obj@active.identDimPlot(obj,label = T)
save(obj, file = 'obj.RData')
###################################### decontX#####################################source('./ks_runDecontX')obj <- ks_runDecontX(seurat_obj = obj,idents = "seurat_clusters")#根据contamination分数删除细胞，至于比例，没有明确的要求，可以按照实际数据处理来设定obj = obj[,obj$estConp < 0.15]


###################################### cell annotation#####################################
featureSets <- list(  "Keratinocytes" = c("KRT5", "KRT10", "KRT14", "KRT15"),   "Fibroblasts" = c("THY1", "COL1A1", "COL11A1"),   "T_cells" = c("CD3D", "CD8A", "CD4","IKZF2", "CCL5"), # PTPRC = CD45   "B_cells" = c("CD19"), # MS41A = CD20, MZB1 = marginal zone B cell specific protein  "APCs" = c("CD14", "CD86", "CD74", "CD163"), # Monocyte lineage (FCGR3A = CD16, FCGR1A = CD64, CD74 = HLA-DR antigens-associated invariant chain)  "Melanocytes" = c("MITF", "SOX10", "MLANA"), # Melanocyte markers  "Endothlial" = c("VWF", "PECAM1", "SELE"), # Endothlial cells (PECAM1 = CD31), SELE = selectin E (found in cytokine stimulated endothelial cells)  "Lymphatic" = c("FLT4", "LYVE1"),  # Lymphatic endothelial (FLT4 = VEGFR-3))  "Muscle" = c("TPM1", "TAGLN", "MYL9"), # TPM = tropomyosin, TAGLN = transgelin (involved crosslinking actin in smooth muscle), MYL = Myosin light chain  "Mast_cells" = c("KIT", "TPSB2", "HPGD"), # KIT = CD117, ENPP3 = CD203c, FCER1A = IgE receptor alpha, TPSB2 is a beta tryptase, which are supposed to be common in mast cells  "HF_surface_markers" = c("ITGB8", "CD200", "SOX9", "LHX2"))
markers <- unlist(featureSets,use.names = F)DotPlot(obj, features = markers,col.min = 0)+coord_flip()table(obj$seurat_clusters)# 0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16   17   18   19 # 2713 2524 1875 1371 1685 1241 1107 1027   64  759  572  480  289  210  212  123   97   78   61   52 
obj <- obj[,!obj@meta.data$seurat_clusters %in% c("8","18","19")]new.cluster.ids <- c("0"="Fibroblasts",                      "1"="Endothlial",                      "2"="Keratinocytes",                      "3"="Myeloid",                      "4"="Keratinocytes",                      "5"="Muscle",                      "6"="Fibroblasts",                      "7"="Muscle",                      "9"="T cells",                      "10"="Keratinocytes",                      "11"="Mast cells",                      "12"="Keratinocytes",                      "13"="Myeloid",                      "14"="Lymphatic",                      "15"="Melanocytes",                      "16"="Keratinocytes",                     "17"="Myeloid")
obj <- RenameIdents(obj, new.cluster.ids)obj$celltype <- obj@active.ident
DimPlot(obj, label = T)+NoLegend()DimPlot(obj, label = T, split.by = 'orig.ident')+NoLegend()
save(obj, file = 'obj.RData')