GSE218208
1.创建Seurat对象
#untar(“GSE218208_RAW.tar”)
rm(list = ls())
a = data.table::fread("GSM6736629_10x-PBMC-1_ds0.1974_CountMatrix.tsv.gz",data.table = F)
a[1:4,1:4]
library(tidyverse)
a$`alias:gene` = str_split(a$`alias:gene`,":",simplify = T)[,1]
#str_split_i(a$`alias:gene`,":",i = 1)
a = distinct(a,`alias:gene`,.keep_all = T) #从数据框a中去除alias:gene列中重复的值，同时保留所有列的信息。
a = column_to_rownames(a,var = "alias:gene") #将数据框a中的alias:gene列的值设置为行名（row names），并将alias:gene列从数据框中移除。
a[1:4,1:4]
library(Seurat)
pbmc <- CreateSeuratObject(counts = a, project = "a", min.cells = 3, min.features = 200)
#使用输入的基因表达矩阵a创建一个新的Seurat对象，并设置项目名称为"a"，同时过滤掉表达在少于3个细胞中的基因，以及过滤掉表达基因数少于200的细胞。
2.质控
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
head(pbmc@meta.data, 3)
VlnPlot(pbmc, features = c("nFeature_RNA","nCount_RNA", "percent.mt"), ncol = 3,pt.size = 0.5)
pbmc = subset(pbmc,nFeature_RNA < 4200 &nCount_RNA < 18000 &percent.mt < 18)
3.降维聚类分群
f = "obj.Rdata"
if(!file.exists(f)){pbmc = pbmc %>% NormalizeData() %>% FindVariableFeatures() %>% ScaleData(features = rownames(.)) %>% RunPCA(pc.genes = pbmc@var.genes) %>%FindNeighbors(dims = 1:15) %>% FindClusters(resolution = 0.5) %>% RunUMAP(dims = 1:15) %>% RunTSNE(dims = 1:15)save(pbmc,file = f)
}
load(f)
ElbowPlot(pbmc)
p1 <- DimPlot(pbmc, reduction = "umap",label = T)+NoLegend();p1
4.makergene
library(dplyr)
f = "markers.Rdata"
if(!file.exists(f)){pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE,min.pct = 0.25)save(pbmc.markers,file = f)
}
load(f)
mks = pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
g = unique(mks$gene)
5.makergene的可视化
DoHeatmap(pbmc, features = g) + NoLegend()+scale_fill_gradientn(colors = c("#2fa1dd", "white", "#f87669"))DotPlot(pbmc, features = g,cols = "RdYlBu") +RotatedAxis()VlnPlot(pbmc, features = g[1:3])FeaturePlot(pbmc, features = g[1:4])
6.注释亚群
手动注释
a = read.delim("../supp/markers.txt",header = F)
gt = split(a[,2],a[,1])DotPlot(pbmc, features = gt,cols = "RdYlBu") +RotatedAxis()
#利用writeLines(paste0(0:11,“,”))，自己手动写，打开一新的text file，将writeLines(paste0(0:11,“,”))的输出写在里边，然后保存在工作目录下，命名为xx.txt
writeLines(paste0(0:11,","))
celltype = read.table("anno.txt",header = F,sep = ",") #自己照着DotPlot图填的
celltype
new.cluster.ids <- celltype$V2
names(new.cluster.ids) <- levels(pbmc)
seu.obj <- RenameIdents(pbmc, new.cluster.ids)
save(seu.obj,file = "seu.obj.Rdata")
p1 <- DimPlot(seu.obj, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
p1
自动注释
library(celldex)
library(SingleR)
ls("package:celldex")
f = "../supp/single_ref/ref_BlueprintEncode.RData"
if(!file.exists(f)){ref <- celldex::BlueprintEncodeData()save(ref,file = f)
}
ref <- get(load(f))
library(BiocParallel)
scRNA = pbmc
test = scRNA@assays$RNA@layers$data
rownames(test) = Features(scRNA)
colnames(test) = Cells(scRNA)
pred.scRNA <- SingleR(test = test, ref = ref,labels = ref$label.main, clusters = scRNA@active.ident)
pred.scRNA$pruned.labels
#查看注释准确性
plotScoreHeatmap(pred.scRNA, clusters=pred.scRNA@rownames, fontsize.row = 9,show_colnames = T)
new.cluster.ids <- pred.scRNA$pruned.labels
names(new.cluster.ids) <- levels(scRNA)
levels(scRNA)
scRNA <- RenameIdents(scRNA,new.cluster.ids)
levels(scRNA)
p2 <- DimPlot(scRNA, reduction = "umap",label = T,pt.size = 0.5) + NoLegend()
p1+p2
可选的celldex包：
a = 1
save(a,file = "a.Rdata")b = load("a.Rdata")b = get(load("a.Rdata")) #load可将a的数值赋值给b
