1.载入包和数据
library(Seurat)
# Load the PBMC dataset 读取数据
pbmc.data <- Read10X(data.dir ="./")
dim(pbmc.data)
#[1] 32738 2700
此时有2700个细胞,32738个features(基因)。
2.创建Seurat对象
pbmc <- CreateSeuratObject(counts = pbmc.data,
project = "pbmc3k", min.cells = 3, min.features = 200)
dim(pbmc)
#[1] 13714 2700
min.cells是过滤features的参数,min.features是过滤细胞的参数。由此可见,导入过程中过滤掉一些features,只剩下13714个基因。
3.数据进行标准化
寻找高变基因之前必须对数据进行标准化,每个细胞的测序深度,或者测序得到的UMI总数是不一样的,因此同一个feature(基因)在不同样品间的表达量是不能直接进行比较的,需要将之标准化。
pbmc <- NormalizeData(pbmc,
normalization.method = "LogNormalize",
scale.factor = 10000)
4.寻找高变基因并可视化
pbmc <- FindVariableFeatures(pbmc,
selection.method = "vst",
nfeatures = 2000)
plot1 <- VariableFeaturePlot(object = pbmc)
plot1
高变基因
在上图中,横坐标为基因在所有细胞中的表达水平(log10对数值),纵坐标为基因在所有细胞中的表达水平的标准差,数值越大,表示该基因在细胞中的表达水平越不稳定。
library(ggplot2)
p <- ggplot(plot1$data,aes(log10(mean),variance.standardized,color = colors))+
geom_point()+theme_bw()
p
高变基因plot
那为什么要挑选高变基因呢,本质其实还是降维(13714个基因降低到2000个),为后续PCA分析,减少计算机运算量。如果后续PCA聚类分析是对所有基因进行,那就没必要寻找高变基因了
那为什么是高变基因而不是恒定表达的基因呢,因为PCA聚类分析中区分不同细胞亚型贡献最大的就是这些高变基因。
网友评论