基础知识
单细胞数据的应用方向
单细胞数据的存放位置
单细胞数据的分析流程
高变基因：方差最大的2000个基因。
marker基因：每个细胞簇中表达显著的基因。
细胞类型注释：根据细胞标志基因进行分类。
具体流程
title: "Seurat"
output: html_document
editor_options:
chunk_output_type: console
knitr::opts_chunk$set(echo = TRUE,warning = F,message = F,fig.width = 10)1.数据和R包准备
代码：https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html
数据：https://s3-us-west-2.amazonaws.com/10x.files/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
代码语言：text复制rm(list = ls())
library(dplyr)
library(Seurat)
library(patchwork)2.读取数据
10X的输入数据是固定的三个文件，在工作目录下新建01_data/，把三个文件放进去。
代码语言：text复制pbmc.data <- Read10X(data.dir = "01_data/")
dim(pbmc.data)
pbmc <- CreateSeuratObject(counts = pbmc.data,
project = "pbmc3k",
min.cells = 3,
min.features = 200)
pbmc
class(pbmc)
# 查看表达矩阵
exp = pbmc@assays$RNA@counts;dim(exp)
exp[30:34,1:4]
#dgCMatrix是稀疏矩阵，点号代表0
range(exp)
boxplot(as.matrix(exp[,1:20]))3.质控
这里是对细胞进行的质控，指标是：
线粒体基因含量不能过高；
nFeature_RNA 不能过高或过低
代码语言：text复制为什么？ nFeature_RNA是每个细胞中检测到的基因数量。nCount_RNA是细胞内检测到的分子总数。nFeature_RNA过低，表示该细胞可能已死/将死或是空液滴。太高的nCount_RNA和/或nFeature_RNA表明“细胞”实际上可以是两个或多个细胞。结合线粒体基因count数除去异常值，即可除去大多数双峰/死细胞/空液滴，因此它们过滤是常见的预处理步骤。 参考自：https://www.biostars.org/p/407036/
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")##符合要求的基因占的百分比；
#f=stringr::str_starts(rownames(pbmc),"MT-")
#rownames(pbmc)[f]
head(pbmc@meta.data)
VlnPlot(pbmc,
features = c("nFeature_RNA",
"nCount_RNA",
"percent.mt"),
ncol = 3,pt.size = 0.5)根据这个三个图，确定了这个数据的过滤标准：
nFeature_RNA在200~2500之间；线粒体基因占比在5%以下。
3.2 三个指标之间的相关性
代码语言：text复制plot1 <- FeatureScatter(pbmc,
feature1 = "nCount_RNA",
feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc,
feature1 = "nCount_RNA",
feature2 = "nFeature_RNA")
plot1 plot23.4 过滤
代码语言：text复制pbmc <- subset(pbmc,
subset = nFeature_RNA > 200 &
nFeature_RNA < 2500 &
percent.mt < 5)
dim(pbmc)4.找高变基因(HVG)
代码语言：text复制pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)
top10 <- head(VariableFeatures(pbmc), 10);top10这里选了2000个，把前十个在图上标记出来。
代码语言：text复制plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1,
points = top10,
repel = TRUE)
plot1 plot25. 标准化和降维
代码语言：text复制all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc[["RNA"]]@scale.data[30:34,1:3]5.1 线性降维PCA
代码语言：text复制pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc))
##只选择了高变化基因分析；
# 查看前5个主成分由哪些feature组成
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
#每个主成分对应基因的热图
DimHeatmap(pbmc, dims = 1:15, cells = 500)
# 应该选多少个主成分进行后续分析
ElbowPlot(pbmc)
# 限速步骤
f = "jc.Rdata"
if(!file.exists(f)){
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
save(pbmc,file = f)
}
load(f)
JackStrawPlot(pbmc, dims = 1:20)#PC1和2
DimPlot(pbmc, reduction = "pca") NoLegend()
# 结合JackStrawPlot和ElbowPlot，挑选10个PC，所以这里dims定义为1:10
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5) #分辨率
#分辨率约高，分的越细致，默认0.5；最大为1，最小为0.1
# 结果聚成几类，用Idents查看
length(levels(Idents(pbmc)))5.2 UMAP 和t-sne
PCA是线性降维，这两个是非线性降维。
代码语言：text复制pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")
#DimPlot(pbmc, reduction = "umap",label=T)6.找marker基因
啥叫marker基因呢。和差异基因里面的上调基因有点类似，某个基因在某一簇细胞里表达量都很高，在其他簇表达量很低，那么这个基因就是这簇细胞的象征。
找全部cluster的maker基因
代码语言：text复制pbmc.markers <- FindAllMarkers(pbmc,
only.pos = TRUE, # 只返回positive基因
min.pct = 0.25) #只计算至少在(两簇细胞总数的)25%的细胞中有表达的基因
#默认参数
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)6.1 比较某个基因在几个cluster之间的表达量
小提琴图
代码语言：text复制VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
VlnPlot(pbmc, features = c("PPBP", "S100A9"))#可以拿count数据画
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)在umap图上标记
代码语言：text复制FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)6.2 marker基因的热图
代码语言：text复制DoHeatmap(pbmc, features = top10$gene) NoLegend()7. 根据marker基因确定细胞
如何去找每种细胞对应的marker基因
注释需要查文献
代码语言：text复制a = read.delim("../supp/markers.txt",header = F)
gt = split(a[,2],a[,1])
DotPlot(pbmc, features = gt,cols = "RdYlBu")
RotatedAxis()new.cluster.ids <- c("Naive CD4 T",
"CD14 Mono",
"Memory CD4 T",
"B",
"CD8 T",
"FCGR3A Mono",
"NK",
"DC",
"Platelet")
names(new.cluster.ids) <- levels(pbmc)
seu.obj <- RenameIdents(pbmc, new.cluster.ids)
p1 <- DimPlot(seu.obj,
reduction = "umap",
label = TRUE,
pt.size = 0.5) NoLegend()
p1
