clean数据经过seeksoultool fast(ffpe默认参数)处理之后,可得到如下matrxi文件:
filtered_feature_bc_matrix
├── barcodes.tsv.gz
├── features.tsv.gz
└── matrix.mtx.gz
该矩阵格式同cellranger等软件处理结果的格式一致,同样适用于下游如Seurat等软件的分析
library(Seurat)
library(tidyverse)
matrix<-Read10X("filtered_feature_bc_matrix")
obj<-CreateSeuratObject(counts=matrix,project = "ffpedemo1")
id_map <- read_delim('filtered_feature_bc_matrix/features.tsv.gz', delim="\t", na="",
col_names = c('Ensembl', 'Symbol')) %>% dplyr::select(1, 2) %>%
mutate(Symbol_uniq=make.unique(Symbol))
rownames(id_map) <- id_map$Symbol_uniq
obj@tools$id_map <- id_map
obj[['percent.mito']] <- PercentageFeatureSet(object = obj, pattern = '^(MT|mt|Mt)-')
obj <- subset(obj, subset = nCount_RNA >= 300)
mt.p <- pnorm(obj$percent.mito, mean = median(obj$percent.mito), sd = mad(obj$percent.mito), lower.tail = FALSE)
mt.lim <- min(obj$percent.mito[which(p.adjust(mt.p, method = "fdr") < 0.05)])
if(mt.lim < 10){mt.lim = 10}
mito_threshold <- mt.lim
if(is.infinite(mito_threshold)){
mito_threshold <- 100
}
obj <- subset(obj, subset = `percent.mito` <= mito_threshold)
obj <- NormalizeData(obj, verbose = FALSE)
obj <- FindVariableFeatures(obj, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
obj <- ScaleData(obj, verbose = FALSE)
obj <- RunPCA(obj, npcs = 15, verbose = FALSE)
obj <- FindNeighbors(obj, dims = 1:15)
resolution <- seq(0.2, 1.4, 0.3)
obj <- FindClusters(obj, resolution = resolution)
resolution_idents <- 'RNA_snn_res.0.8'
Idents(obj) <- resolution_idents
obj[['seurat_clusters']] <- obj[[resolution_idents]]
obj <- RunUMAP(obj, reduction = "pca", dims = 1:15)
obj <- RunTSNE(obj, reduction = "pca", dims = 1:15, check_duplicates = FALSE)
saveRDS(obj, file="ffpedemo1.rds")
ffpe<-readRDS("ffpedemo1.rds")
options(repr.plot.height=6, repr.plot.width=7)
DimPlot(ffpe, reduction = "umap", label = TRUE)
options(warn = -1)
options(repr.plot.height=4, repr.plot.width=8)
FeaturePlot(ffpe, features = c("nFeature_RNA","nCount_RNA"))
VlnPlot(ffpe, features = c("nFeature_RNA","nCount_RNA"))
VlnPlot(ffpe,features = c("percent.mito"))
## anno
options(repr.plot.height=6, repr.plot.width=20)
DotPlot(ffpe,features = c(
"SFTPB","MUC1", #Alveolar cell type 2
"EPCAM", "KRT18", "KRT19", #Epithelial
"S100A8","S100A9", # Neutrophil
"CD3E","CD3D", # Tcell
"PLP1","CTNNA3","IL1RAPL1",#Oligodendrocytes
"OLIG1","LHFPL3" ,"PDGFRA", "COL9A1", #"PTPRZ1",# OPC
"SLC1A2","GFAP", "SLCO1C1", "AQP4", # Astrocytes
"PTPRC", "P2RY12", "ITGAM", #Microglia
"NEUROD2", "RBFOX1", "SLC17A7",#Ex-Neuron
"CLDN5", "RGS5" ,#Endothelial
"GAD1","GAD2","PDE4DIP","SLC32A1","DLX1",#InterNeurons
"PTN", "COL1A2" #Pericyte
)) + theme(axis.text.x = element_text(angle=60,vjust = 0.6))
annoffpe1 <- RenameIdents(ffpe,
"0"="Oligodendrocytes",
"1"="Epithelial",
"2"="Epithelial",
"3"="Epithelial",
"4"="Microglia",
"5"="OPC",
"6"="Oligodendrocytes",
"7"="Ex-Neuron",
"8"="Pericyte",
"9"="Microglia",
"10"="Neutrophil",
"11"="Microglia",
"12"="Tcell",
"13"="Astrocytes",
"14"="Epithelial",
"15"="InterNeurons",
"16"="Endothelial"
)
annoffpe1$cell_type<-Idents(annoffpe1)
options(repr.plot.height=6, repr.plot.width=10)
DimPlot(annoffpe1, reduction = "umap", label = TRUE)