单细胞数据分析2（练习GSE218208)

原创

Erics blog

发布于 2023-10-15 16:02:13

4140

发布于 2023-10-15 16:02:13

文章被收录于专栏：R语言数据分析R语言数据分析

title: "GSE218208"

output: html_document

editor_options:

chunk_output_type: console

knitr::opts_chunk$set(echo = TRUE,warning = F,message = F,fig.width = 10)

1.创建Seurat对象

#https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10039079/
#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 = column_to_rownames(a,var = "alias:gene")
a[1:4,1:4]
#从GEO下载的数据需要自己处理后读取
library(Seurat)
pbmc <- CreateSeuratObject(counts = a, 
                           project = "a", ##自己命名
                           min.cells = 3, 
                           min.features = 200)
exp = pbmc[["RNA"]]@counts;dim(exp)
exp[1:4,1:4]

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 %>% 
  Seurat::NormalizeData() %>%  
  FindVariableFeatures() %>%  
  ScaleData() %>%  
  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, logfc.threshold = 0.25)
  save(pbmc.markers,file = f)
}
load(f)
mks = pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
mks
g = unique(mks$gene)
g

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)

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,","))
celltype = read.table("anno.txt",header = F,sep = ",") #自己照着DotPlot图填的
#R语言可以直接新建txt文档
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