背景
用R做单细胞RNA-seq的t-SNE图。
下图来自文章PMID: 28726821：
一般可能会卡在两个问题上：R包安装和数据集的预处理。
本代码带你顺利解决这两个问题，画出t-SNE图，把细胞分类和候选marker基因输出到文件。
应用场景
场景一：单细胞测序的上百个细胞，识别亚群，找到marker基因；
场景二：上百个肿瘤样本，识别亚型，找到marker基因。像例图那样，用来划分甲基化亚型。
安装R包 Seurat
# loading Seurat
if ( ! require("Seurat")){
install.packages('Seurat')
library(Seurat)
} else{
library(Seurat)
}
# loading plyr
if (! require(plyr)){
install.packages("plyr")
library(plyr)
} else{
library(plyr)
}
注: Windows系统用户安装会很顺利。
MacOS/Linux系统用户要有Root权限，或许会因为hdf5r安装失败而无法顺利安装，不要图省事，按报错提示，用 brew install hdf5。
因为根据hdf5的源代码，它不支持最新的hdf5 1.10.2，最高只支持hdf5 1.10.0，所以请在root权限下根据源码安装。
wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.0/src/hdf5-1.10.0.tar.gz
tar xf hdf5-1.10.0.tar.gz
cd hdf5-1.10.0
./configure --prefix=/usr/local
make && sudo make install
数据预处理
请先解压当前目录的下的20250408.zip, 里面存放着测试数据集
输入数据是基于UMI的表达矩阵而不是count数据，其中数据集中的行表示基因，列表示细胞。
来自文章PMID: 28045081
读取数据
molecules <- read.delim("./tung/molecules.txt",header = TRUE,
row.names = 1)
dim(molecules)
#查看前10个基因在前3个细胞中的表达量
molecules[1:10,1:3]
molecules前几行如下：
创建Seurat对象
数据初始化，并初步筛选
library(Seurat)
obj <- CreateSeuratObject(counts = molecules,
min.cells = 3, #筛选至少在3个细胞中表达的基因
min.features = 200) #筛选至少有200个基因表达的细胞
增加Seurat的元信息(metadata)
元信息存放在 [email protected]中，可以是细胞的分组信息，所做处理等，需要自己提供。
此处包含3个人，各96个细胞，各3次重复。
ann <- read.delim("annotation.txt",header = TRUE)
rownames(ann) <- ann$sample_id
head(ann)
summary(ann)
ann前几行如下：
检查样品信息跟表达量矩阵的样品是否一致
if (all(rownames(ann) %in% [email protected])){
obj <- AddMetaData(obj, ann)
} else{
warning("row names should be same to the items in [email protected]")
}
过滤掉低质量细胞
根据具体情况来确定标准，进行QC，过滤掉低质量的细胞，此处用 VlnPlot探索可能的离群值
VlnPlot(obj, features = c("nCount_RNA","nFeature_RNA"), group.by = "orig.ident")
一般我们会过滤掉过低表达基因的细胞，上图中，你可以认为nCount_RNA低于 25000、nFeature_RNA低于 6000的细胞是低质量细胞。把这两个值写到下面过滤细胞的参数中：
obj <- subset(obj, subset = nCount_RNA > 25000 & nFeature_RNA > 6000)
此处所用的数据实验深度足够，10X技术单细胞数据标准会不同。
注意：笔者用的‘Seurat’ version 5.2.1，新版语法和旧版本不同。
标准化
标准化这一步和混池RNA-seq类似，一般不会只使用raw count，而会用log标准化等方式降低其他干扰因素的影响。
obj <- NormalizeData(
object = obj,
normalization.method = "LogNormalize",
scale.factor = 10000
)
寻找细胞间变化明显的基因
直接使用全部基因进行后续分析对配置要求很高，而是选择那些在细胞中具有明显变化的基因。在降低资源消耗的同时还保证了结果的可靠性。
obj <- FindVariableFeatures(
object = obj,
selection.method = "vst", # 默认方法，替代旧版参数
nfeatures = 2000, # 指定高变基因数量（默认2000）
mean.cutoff = c(0.0125, 3),# 替代旧版 x.low.cutoff 和 x.high.cutoff
dispersion.cutoff = 0.5 # 替代旧版 y.cutoff
)
处理混淆因子
所谓的混淆因子(confounder factor) 指的就是目标因素（如实验处理）以外的因素对结果的影响， 比如说批次效应。这些因素掩盖了你真正想研究的目标，所以要尽可能排除。
obj <- ScaleData(
object = obj,
vars.to.regress = c("nCount_RNA")
)
线性降维（PCA）
后续会对细胞进行聚类，直接使用所有特征（即基因）会出现数据挖掘中的“维度爆炸”，即维度越高反而结果越差。因此我们需要用PCA进行线性降维，使用部分主成分代替所有基因进行后续的t-SNE和细胞聚类分析。
根据PCElbowPlot确定主成分个数，曲线变得平缓的横坐标。
这个主成分个数将用在下一部分的 pca_to_use =参数中
obj <- RunPCA(
object = obj,
pc.genes = [email protected],
do.print = FALSE
)
ElbowPlot(object = obj)
从图中可以看出：横坐标PC为14时，曲线趋近平缓，因此，此处定为14。
细胞聚类
细胞聚类是单细胞分析的关键一步，毕竟大部分研究是想找到新的细胞分群、以前没有发现的细胞类型
# 构建 SNN 图（关键步骤）
pca_to_use <- 14 # 上一步的数目
obj <- FindNeighbors(obj, dims = 1:pca_to_use, reduction = "pca") # 默认生成 RNA_snn 图
# 新版语法（正确方式）
obj <- FindClusters(
object = obj,
graph.name = "RNA_snn", # 默认图名（来自 FindNeighbors）
resolution = 1.0 # 分辨率控制分群数
)
# resolution：决定分群数目，值越大分群越细（如 0.8 分群数多于 0.5）。
# graph.name：若使用自定义图名（如 integrated_snn），需显式指定
# 统计各簇细胞数
table(obj$seurat_clusters)
# 查看元数据中的分群信息
head([email protected])
t-SNE非线性降维
obj <- RunTSNE(
object = obj,
dims.use = 1:pca_to_use,
do.fast = TRUE
)
开始画图
按照细胞分群(group.by="ident")进行着色
设置 pt_size修改点的大小，修改 mycol增加自定义颜色
#定义足够多的颜色，后面从这里选颜色
mycol <- c("#223D6C","#D20A13","#FFD121","#088247","#11AA4D","#58CDD9","#7A142C","#5D90BA","#431A3D","#91612D","#6E568C","#E0367A","#D8D155","#64495D","#7CC767")
t-SNE可视化
pt_size = 1
df <- obj@[email protected]
df <- cbind(df, as.data.frame(obj@ident))
colnames(df) <- c("x", "y", "ident")
p <- ggplot(df, aes(x=x, y=y)) + geom_point(aes(colour=ident),size=pt_size) +
labs(x="t-SNE 1", y ="t-SNE 2")
画椭圆线
用于辅助区分不同细胞类群
obs <- p$data[,c("x","y","ident")]
ellipse_pro <- 0.98
theta <- c(seq(-pi, pi, length = 50), seq(pi, -pi, length = 50))
circle <- cbind(cos(theta), sin(theta))
ell <- ddply(obs, 'ident', function(x) {
if(nrow(x) <= 2) {
return(NULL)
}
sigma <- var(cbind(x$x, x$y))
mu <- c(mean(x$x), mean(x$y))
ed <- sqrt(qchisq(ellipse_pro, df = 2))
data.frame(sweep(circle %*% chol(sigma) * ed, 2, mu, FUN = '+'))
})
names(ell)[2:3] <- c('x', 'y')
ell <- ddply(ell, .(ident) , function(x) x[chull(x$x, x$y), ])
p1 <- p + geom_polygon(data = ell, aes(group = ident),
colour = "black",
alpha = 1,fill = NA,
linetype="dashed")
# linetype可选类型:
# blank, solid, dashed, dotted, dotdash, longdash, twodash
发表于 2025-4-8 14:46:14
|
查看: 116|
回复: 2
|Archiver|手机版|小黑屋|bioinfoer
( 萌ICP备20244422号 )
