|
我們介紹的是scATAC聯(lián)合scRNA數(shù)據(jù)分析，那么scRNA必不可少，我們選擇的數(shù)據(jù)相同的分組也測了scRNA，所以還是回顧一下，跑一下scRNA流程，做好細胞分群和注釋吧，一方面注釋好的scRNA可用于scATAC數(shù)據(jù)注釋，另一方面兩者可以聯(lián)合分析做一些內(nèi)容。話不多說，都很熟悉，直接上代碼：上游分析不再多說！
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")#根據(jù)contamination分數(shù)刪除細胞，至于比例，沒有明確的要求，可以按照實際數(shù)據(jù)處理來設(shè)定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')
|
