library(Seurat)
library(dplyr)
library(magrittr)
library(patchwork)
library(hdf5r)
library(tidyverse)
model.data <- Read10X(data.dir = "/home/jaa20/bak/jaa20/model/Model/outs/count/filtered_feature_bc_matrix")
XHP10W.data <- Read10X(data.dir = "/home/jaa20/bak/jaa20/XHP10W/XHP10W/outs/count/filtered_feature_bc_matrix")
创建Seurat对象
使用CreateSeuratObject函数创建Seurat对象，这里要求细胞中基因数目不能小于200，且基因至少在三个细胞中有表达，否则过滤掉
model_f <- CreateSeuratObject(counts = model.data, project = "modelf", min.cells = 3, min.features = 200)
XHP10W_f <- CreateSeuratObject(counts = XHP10W.data, project = "XHP10Wf", min.cells = 3, min.features = 200)
Seurat对象就是一个S4类，里面装着单细胞数据集，如count矩阵、细胞矩阵、聚类信息都存储在这个容器中。
model_f[["percent.mt"]] <- PercentageFeatureSet(model_f, pattern = "^MT-")
VlnPlot(model_f, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
XHP10W_f[["percent.mt"]] <- PercentageFeatureSet(XHP10W_f, pattern = "^MT-")
VlnPlot(XHP10W_f, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
我们过滤具有超过2500或少于200个基因计数的细胞，同时过滤掉线粒体比例超过5%的细胞。
model_f1 <- subset(model_f, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
XHP10W_f1 <- subset(XHP10W_f, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
下一步是标准化数据。默认情况下，我们使用“lognormalize”全局缩放的归一化方法，通过总表达值对每个细胞的基因表达值归一化，并将其乘以缩放因子(默认为10,000)，最后对结果进行对数变换
model_f1 <- NormalizeData(model_f1, normalization.method = "LogNormalize", scale.factor = 10000)
XHP10W_f1 <- NormalizeData(XHP10W_f1, normalization.method = "LogNormalize", scale.factor = 10000)
计算数据集中表现出高细胞间差异的基因子集(即，它们在一些细胞中高表达，而在另一些细胞中低表达)。使用均值与方差之间的关系，来挑选高变基因，在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。默认情况下，每个数据集选择2000个高变异基因用于下游分析。
model_f1 <- FindVariableFeatures(model_f1, selection.method = "vst", nfeatures = 2000)
XHP10W_f1 <- FindVariableFeatures(XHP10W_f1, selection.method = "vst", nfeatures = 2000)
model_f1 <- ScaleData(model_f1)
XHP10W_f1 <- ScaleData(XHP10W_f1)
VlnPlot(model_f1,
features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB"),
cols =rainbow(col.num),
pt.size = 0.1,
ncol = 4) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
VlnPlot(XHP10W_f1,
features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB"),
cols =rainbow(col.num),
pt.size = 0.1,
ncol = 4) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
上一步找到的高变基因，常常会包含一些细胞周期相关基因。它们会导致细胞聚类发生一定的偏移，即相同类型的细胞在聚类时会因为细胞周期的不同而分开。
细胞周期评分
g2m_genes_model = cc.geness.genes
s_genes_model = CaseMatch(search = s_genes_model, match = rownames(model_f1))
model_f1 <- CellCycleScoring(object=model_f1, g2m.features=g2m_genes_model, s.features=s_genes_model)
查看细胞周期基因对细胞聚类的影响
model_f1a <- RunPCA(model_f1, features = c(s_genes_model, g2m_genes_model))
p_model <- DimPlot(model_f1a, reduction = "pca", group.by = "Phase")
ggsave("p_model_pca.png", p_model, width = 8, height = 6)
如果需要消除细胞周期的影响
model_f11 <- ScaleData(model_f1, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(model_f1))
g2m_genes_model = cc.genes$g2m.genes
g2m_genes_model = CaseMatch(search = g2m_genes_model, match = rownames(model_f11))
s_genes_model = cc.genes$s.genes
s_genes_model = CaseMatch(search = s_genes_model, match = rownames(model_f11))
model_f11 <- CellCycleScoring(object=model_f1, g2m.features=g2m_genes_model, s.features=s_genes_model)
##查看细胞周期基因对细胞聚类的影响
model_f11a <- RunPCA(model_f11, features = c(s_genes_model, g2m_genes_model))
p_model1 <- DimPlot(model_f11a, reduction = "pca", group.by = "Phase")
ggsave("p_model1_pca.png", p_model1, width = 8, height = 6)
g2m_genes_XHP10W = cc.geness.genes
s_genes_XHP10W = CaseMatch(search = s_genes_XHP10W, match = rownames(XHP10W_f1))
XHP10W_f1 <- CellCycleScoring(object=XHP10W_f1, g2m.features=g2m_genes_XHP10W, s.features=s_genes_XHP10W)
查看细胞周期基因对细胞聚类的影响
XHP10W_f1a <- RunPCA(XHP10W_f1, features = c(s_genes_XHP10W, g2m_genes_XHP10W))
p_XHP10W <- DimPlot(XHP10W_f1a, reduction = "pca", group.by = "Phase")
ggsave("p_XHP10W_pca.png", p_XHP10W, width = 8, height = 6)
