幸运快三官方首页

胃癌单细胞数据集GSE163558复现(八):拟时序分析
发布日期:2025-03-07 16:41    点击次数:56

图片

前言

Hello小伙伴们大家好,我是生信技能树的小学徒”我才不吃蛋黄“。今天是胃癌单细胞数据集GSE163558复现系列第八期。第七期我们使用AddModuleScore_UCell函数计算细胞的增殖和迁移评分。本期,我们将使用monocle2进行拟时序分析。

1.背景介绍

系列推文前七期,我们共同学习了单细胞测序的基础分析,从第八期开始,我们将陆续学习拟时序分析(Pseudo-Temporal Analysis),“Copykat”分析及细胞通讯等高级分析。

疾病是一个动态变化的过程,在疾病发生发展的过程中,细胞状态随时间不断变化。现有的单细胞测序技术提供的是某一时刻的“快照”,然而通过拟时序分析,我们可以推断出细胞的动态变化过程。

拟时序分析(Pseudo-Temporal Analysis)的定义如下:

拟时序分析(Pseudo-Temporal Analysis)是一种分析方法,通常用于研究那些在时间上不连续的数据集。它试图通过模拟时间序列的方式来分析数据,即使数据本身并不是按照时间顺序收集的。这种分析方法在多个领域都有应用,比如在生物信息学中,它可以用于分析基因表达数据,即使这些数据不是在连续的时间点上收集的。拟时序分析的关键在于构建一个虚拟的时间轴,然后在这个时间轴上对数据进行排序和分析。这通常涉及到以下几个步骤:数据预处理:清洗数据,确保数据的质量和一致性。特征选择:识别和选择对分析有用的特征或变量。时间轴构建:根据数据的特性构建一个虚拟的时间轴,这可能基于生物学过程、实验设计或其他逻辑。排序:将数据按照虚拟时间轴进行排序。分析:使用时间序列分析的方法来分析排序后的数据,比如趋势分析、周期性分析等。验证:通过交叉验证或其他方法来验证分析结果的可靠性。*

单细胞拟时序分析的主要目的包括:

细胞发展轨迹推断:通过分析细胞状态的变化,推断细胞是如何从一种状态发展或分化到另一种状态。细胞状态排序:将细胞按照它们在假定的发展轨迹上的位置进行排序,即使这些细胞并非在连续的时间点上被观测。生物学过程理解:识别细胞在特定生物学过程中的动态变化,如细胞分化、发育、疾病进展等。关键基因和调控网络识别:发现在细胞状态转变中起关键作用的基因和调控网络。细胞异质性探究:揭示细胞群体内部的异质性,以及不同细胞亚群之间的关系。细胞命运预测:基于当前的细胞状态和拟时序轨迹,预测细胞可能的未来发展和分化命运。动态生物学过程建模:构建数学模型来模拟细胞状态的动态变化,为生物学假设提供定量描述。疾病研究和治疗靶点发现:在疾病研究中,拟时序分析可以帮助理解病理状态下细胞的变化,发现潜在的治疗靶点。发育过程的细胞谱系构建:在发育生物学中,拟时序分析有助于构建细胞谱系树,理解细胞的起源和命运。实验设计指导:通过理解细胞状态的变化,可以指导未来的实验设计,例如选择哪些细胞状态进行进一步的实验验证。

Monocle是目前最常用的拟时序分析算法,其核心技术是一种机器学习算法——反向图形嵌入(Reversed Graph Embedding)进行降维分析。Monocle2使用DDRtree进行降维分析,而Monocle3则使用UMAP降维,可视化细胞转录特征相似性关系,从而描述细胞状态过渡轨迹。细胞状态过渡轨迹是“根”到“叶”的树状结构,细胞从根到达分支后,会选择一个分支往远处移动,最后到达叶子,每个细胞的psudotime值是它从叶返回根的距离。

monocle的分析流程与seurat类似,准备输入文件,创建monocle对象cds(new cell data set),对数据预处理,筛选特定基因,降维,计算psudotime值,最后可视化。本期使用monocle2进行拟时序分析。

2.数据分析2.1 导入数据

清除系统环境变量,设置工作目录,加载R包,选择做拟时序的亚群,读取恶性上皮细胞RDS文件:

rm(list=ls())getwd()setwd('')library(tidyverse)library(tinyarray)library(data.table) library(Seurat)library(ggplot2)library(clustree)library(cowplot)library(dplyr)library(monocle)sce.all = readRDS('malignant.rds')sce = sce.alltable(sce$celltype)Idents(sce) = sce$celltype

如果用自己电脑,细胞量太大,可以每个细胞亚群抽样 :

allCells=names(Idents(sce))allType = levels(Idents(sce))# choose_Cells = unlist(lapply(allType, function(x){#   cgCells = allCells[Idents(sce)== x ]#   cg=sample(cgCells,10)#   cg# }))# cg_sce = sce[, allCells %in% choose_Cells]cg_sce = scetable(Idents(cg_sce))
2.2 数据准备

monocle构建CDS需要3个矩阵:

expr.matrixphenodata(pd)featuredata(fd)

expr.matrix:

Mono_tj<-cg_sceMono_matrix<-as(as.matrix(GetAssayData(Mono_tj,slot = "counts")), 'sparseMatrix')

构建featuredata:

一般需要两个col,一个是gene_id,一个是gene_short_name,row对应counts的rownames:

feature_ann<-data.frame(gene_id=rownames(Mono_matrix),gene_short_name=rownames(Mono_matrix))rownames(feature_ann)<-rownames(Mono_matrix)Mono_fd<-new("AnnotatedDataFrame", data = feature_ann)

构建phenodata:

Seurat object中的@meta.data一般会存放表型相关的信息如cluster、sample的来源、group等,所以选择将metadata转换为phenodata:

sample_ann<-Mono_tj@meta.datarownames(sample_ann)<-colnames(Mono_matrix)Mono_pd<-new("AnnotatedDataFrame", data =sample_ann)
2.3 monocle分析流程

构建new cell data set:

Mono.cds<-newCellDataSet(Mono_matrix,phenoData =Mono_pd,featureData =Mono_fd,expressionFamily=negbinomial.size())

newCellDataSet函数中,expressionFamily参数用于指定表达矩阵的数据类型,有几个选项可以选择:稀疏矩阵用negbinomial.size();FPKM值用tobit();logFPKM值用gaussianff()。

查看phenodata、featuredata:

head(pData(Mono.cds))head(fData(Mono.cds))

Mono.cds相当于Seurat V5构建的Seurat对象,随后需要对数据预处理(估计尺度因子、估计离散度):

Mono.cds <- estimateSizeFactors(Mono.cds)Mono.cds <- estimateDispersions(Mono.cds)

筛选基因,这里可以根据自己的需要筛选特定的基因:

disp_table <- dispersionTable(Mono.cds)unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)Mono.cds <- setOrderingFilter(Mono.cds, unsup_clustering_genes$gene_id)

用DDRtree进行降维分析:

Mono.cds <- reduceDimension(  Mono.cds,  max_components = 2,  method = 'DDRTree')

计算psudotime值:

Mono.cds <- orderCells(Mono.cds)head(pData(Mono.cds))

图片

2.4 monocle可视化

展示State轨迹分布图:

plot_cell_trajectory(Mono.cds,cell_size = 1)

图片

展示Cluster/Pseudotime轨迹分布图:

p1 = plot_cell_trajectory(Mono.cds,color_by="celltype", size=1,show_backbone=TRUE)p2 = plot_cell_trajectory(Mono.cds,color_by="Pseudotime", size=1,show_backbone=TRUE) p1+p2

图片

分面显示:

plot_cell_trajectory(Mono.cds, color_by = "celltype") + facet_wrap("~sample", nrow = 1)

图片

树形图:

plot_complex_cell_trajectory(Mono.cds,x=1,y=2,color_by="celltype")+  scale_color_manual(values =mycolors)+  theme(legend.title = element_blank())

图片

沿时间轴的细胞密度图:

library(ggpubr)df <- pData(Mono.cds)view(df)

图片

ggplot(df,aes(Pseudotime, colour = celltype, fill=celltype))+        geom_density(bw=0.5,size=1,alpha =0.5)+theme_classic2()

图片

Monocle基因可视化:

head(unsup_clustering_genes)s.genes <- c("NOC2L","PLEKHN1","HES4","ISG15")p1 <- plot_genes_jitter(Mono.cds[s.genes,],                        grouping = "State", color_by = "State")p2 <- plot_genes_violin(Mono.cds[s.genes,],                        grouping = "State", color_by = "State")p3 <- plot_genes_in_pseudotime(Mono.cds[s.genes,], color_by = "State")

图片

图片

图片

拟时相关基因聚类热图:

#高变基因disp.genes <- subset(disp_table, mean_expression >= 0.5&dispersion_empirical >= 1*dispersion_fit)disp.genes <- as.character(disp.genes$gene_id)diff_test <- differentialGeneTest(Mono.cds[disp.genes,], cores = 4,                                   fullModelFormulaStr = "~sm.ns(Pseudotime)")sig_gene_names <- row.names(subset(diff_test, qval < 1e-50))plot_pseudotime_heatmap(Mono.cds[sig_gene_names,], num_clusters=4,                             show_rownames=T, return_heatmap=T)

图片

2.5 BEAM分析

由于细胞基因表达模式不同,单细胞轨迹中通常包括分支。通过BEAM(Branched expression analysis modeling)分析,我们可找到以依赖于分支的方式调控的基因。

disp_table <- dispersionTable(Mono.cds)disp.genes <- subset(disp_table, mean_expression >= 0.5&dispersion_empirical >= 1*dispersion_fit)disp.genes <- as.character(disp.genes$gene_id)mycds_sub <- Mono.cds[disp.genes,]plot_cell_trajectory(mycds_sub, color_by = "State")beam_res <- BEAM(mycds_sub, branch_point = 1, progenitor_method = "duplicate")beam_res <- beam_res[order(beam_res$qval),]beam_res <- beam_res[,c("gene_short_name", "pval", "qval")]mycds_sub_beam <- mycds_sub[row.names(subset(beam_res, qval < 1e-4)),]plot_genes_branched_heatmap(mycds_sub_beam,  branch_point = 1, num_clusters = 3, show_rownames = T)

图片

结语

本期,我们进入了单细胞高级分析系列,使用monocle2进行了拟时序分析。下一期,根据课程安排,我们将暂时回到Bulk转录组,利用TCGA-STAD数据进行生存分析。干货满满,欢迎大家持续追更,谢谢!

图片

本站仅提供存储服务,所有内容均由用户发布,如发现有害或侵权内容,请点击举报。

友情链接: