当前位置：首页 > news >正文

31-植物单细胞RNA-seq分析教程5-2025年版

news 2026/3/27 2:16:45

本教程基于银白杨顶芽单细胞测序数据，聚焦 CYC 和 CDK 基因家族参与叶片发育的分子机制，适配有少量生信基础（了解 Linux 基本操作、R 语言入门）的读者。教程共4节，遵循 “软件准备→基础分析→细胞注释→高级挖掘” 的分析流程，代码均来自实际研究，注释结合发表文章核心结果，确保 “代码可复现、结果可解读”。

第4节，拟时序分析

####拟时序分析是模拟各个组织之间的发育过程，这里需要根据差异基因进行分析，不同差异基因、不同的软件、不同的参数得到的轨迹是天差地别的，如使用seurat选择的高变基因，使用clusters差异表达基因，使用monocle选择的高变基因。另外细胞组织类型的选择，也有很大的影响。我建议大家都把主流差异基因计算的轨迹都算一次，算出20个左右，然后根据真实的发育规律判断哪一个轨迹适合放到文章里，如叶片一定从叶原基起源，叶原基分布在茎尖分生组织SAM中，所以SAM一定是轨迹开端，需要根据实际情况而选择。#######

#########我文章中选择的是monocle的差异基因，含有未鉴定的组织类群##############
conda activate R
R
library(Seurat);library(tidyverse);library(patchwork);library(dplyr);library(ggplot2);library(monocle)
load("all.scRNA1.Rdata")
#注释细胞类型
new.cluster.ids<- c("0","1","SMC","MC","TRI","EC","IMC","VC","8","MC","PC","MC","EC","IMC","SMC","EC","VC")
names(new.cluster.ids) <- levels(scRNA1)
scRNA1<- RenameIdents(scRNA1,new.cluster.ids)
scRNA1[['cell_type']]<- scRNA1@active.ident
scRNA1_subset <- subset(scRNA1, idents = c("SMC","MC","TRI","EC","IMC","VC","PC","0","1","8"))
scRNA1_subset <- subset(x = scRNA1_subset, cells = sample(x = 1:length(Idents(scRNA1_subset)), size = round(0.3 * length(Idents(scRNA1_subset))), replace = FALSE))
scRNA1_subset#30167 features across 8703 samples
table(Idents(scRNA1_subset))
#   0    1  SMC   MC  TRI   EC  IMC   VC    8   PC
#1240 1067 1084 1565  646  960  666  625  478  372
#构建对象
cds<- as.CellDataSet(scRNA1_subset)
cds<- estimateSizeFactors(cds)
cds<- estimateDispersions(cds)disp_table <- dispersionTable(cds)
disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical>=1 * dispersion_fit)$gene_id
cds <- setOrderingFilter(cds, disp.genes)
cds <- reduceDimension(cds, max_components = 2, method = 'DDRTree')
cds <- orderCells(cds,reverse=T)
length(disp.genes)#1277
diff<- differentialGeneTest(cds[disp.genes,],fullModelFormulaStr = "~cell_type")
head(diff)
deg <- subset(diff, qval < 0.01)
deg <- deg[order(deg$qval,decreasing=F),]#排序
head(deg)
write.table(deg,file="train.monocle.DEG.xls",col.names=T,row.names=F,sep="\t",quote=F)
ordergene <- rownames(deg)
cds<- setOrderingFilter(cds, ordergene)
pdf("ordering_genes_cds_2.pdf",height=5,width=10)
plot_ordering_genes(cds)
dev.off()
#1，以pseudotime值上色
pdf("train.monocle.pseudotime.pdf",width = 7,height = 7)
plot_cell_trajectory(cds,color_by="Pseudotime", cell_size=0.5,show_backbone=TRUE) 
dev.off()
#2，以细胞类型上色
pdf("train.monocle.celltype.pdf",width = 7,height = 7)
plot_cell_trajectory(cds,color_by="cell_type", cell_size=0.5,show_backbone=TRUE)
dev.off()
#2.1，可以看每个cluster的分布情况
pdf("train.monocle.celltype_everyone.pdf",width = 12,height = 10)
plot_cell_trajectory(cds,color_by="cell_type", cell_size=0.5,show_backbone=TRUE) + facet_wrap("~cell_type", nrow = 2, ncol=5)
dev.off()
#3，以细胞状态上色
pdf("train.monocle.state.pdf",width = 7,height = 7)
plot_cell_trajectory(cds, color_by = "State",cell_size=0.5,show_backbone=TRUE)
dev.off()
#4，以seurat分群排序细胞
pdf("train.monocle.seurat.clusters.pdf",width = 7,height = 7)
plot_cell_trajectory(cds, color_by = "seurat_clusters",cell_size=0.5,show_backbone=TRUE)
dev.off()
#5，画树状图
pdf("train.monocle.tree.pdf",width = 7,height = 7)
plot_complex_cell_trajectory(cds, x = 1, y = 2,color_by = "cell_type",cell_size=0.5)
dev.off()
#6，细胞密度图
library("ggpubr")
pdf("train.monocle.density.pdf",width = 7,height = 7)
df<- pData(cds)#pData(cds)取出的是cds对象中cds@phenoData@data的内容
ggplot(df, aes(Pseudotime, colour = cell_type, fill=cell_type)) +geom_density(bw=0.5,size=1,alpha = 0.5)+theme_classic2()
dev.off()
####指定基因的可视化
keygenes<- head(ordergene,4)
cds_subset <- cds[keygenes,]
p1 <- plot_genes_in_pseudotime(cds_subset, color_by = "State")
p2 <- plot_genes_in_pseudotime(cds_subset, color_by = "cell_type")
p3 <- plot_genes_in_pseudotime(cds_subset, color_by = "Pseudotime")
plotc <- p1|p2|p3
ggsave("Genes_pseudotimeplot.pdf", plot = plotc, width = 16, height = 8)#######以下是一些特殊基因单独拎出来的分析##################
#指定基因CYCD1;4, CYCD3;3, CYCD3;5, CYCB1;1, CDKA;1, CDKB2;2
s.genes <- c("evm.TU.Poalb09G007780.1","evm.TU.Poalb05G011570.1","evm.TU.Poalb14G002050.1", "evm.TU.Poalb16G006810.1","evm.TU.Poalb04G011340.1","evm.TU.Poalb05G021370.1")
p1 <- plot_genes_jitter(cds[s.genes,], grouping = "State", color_by = "State")
p2 <- plot_genes_violin(cds[s.genes,], grouping = "State", color_by = "State")
p3 <- plot_genes_in_pseudotime(cds[s.genes,], color_by = "State")
plotc <- p1|p2|p3
ggsave("Genes_Jitterplot.pdf", plot = plotc, width = 16, height = 8)
#拟时序展示单个基因的表达量CYCD3;3, CYCD3;5
colnames(pData(cds))
pData(cds)$evm.TU.Poalb05G011570.1<- log2(exprs(cds)['evm.TU.Poalb05G011570.1',]+1)
p1<- plot_cell_trajectory(cds, color_by="evm.TU.Poalb05G011570.1")
pData(cds)$evm.TU.Poalb14G002050.1<- log2(exprs(cds)['evm.TU.Poalb14G002050.1',]+1)
p2<- plot_cell_trajectory(cds, color_by="evm.TU.Poalb14G002050.1")
plotc<- p1+p2
ggsave("Genes_pseudotimeplot2.pdf", plot = plotc, width = 8, height = 5)
####寻找拟时序相关的基因
ordergene<- ordergene[grep("\\.1$", ordergene)]
ordergene <- ordergene[!grepl("\\.1\\.1$", ordergene)]
ordergene <- ordergene[grep("^evm.TU", ordergene)]
Time_diff <- differentialGeneTest(cds[ordergene,], fullModelFormulaStr = "~sm.ns(Pseudotime)")
write.csv(Time_diff,"Time_diff_all.csv",row.names=F)
Time_genes <- Time_diff %>% pull(gene_short_name) %>% as.character()
p<- plot_pseudotime_heatmap(cds[Time_genes,], num_clusters=4, show_rownames=F, return_heatmap=T)
ggsave("Time_heatmapAll.pdf",p, width=10, height=10)
#以下是以拟时差异基因热图绘制。拟时差异基因热图绘制（提取了前100个）
Time_genes <- top_n(Time_diff, n = 100, desc(qval)) %>% pull(gene_short_name) %>% as.character()
p <- plot_pseudotime_heatmap(cds[Time_genes,], num_clusters=2, show_rownames=F, return_heatmap=T)#聚类两列，不展示名字。
ggsave("Time_heatmapTop100.pdf", p, width = 8, height = 8)
hp.genes <- p$tree_row$labels[p$tree_row$order]
Time_diff_sig <- Time_diff[hp.genes, c("gene_short_name", "pval", "qval")]
write.csv(Time_diff_sig, "Time_diff_sig.csv", row.names = F)#保存，这个图没什么用
####保存单个变量，拟时序画图都是用CDS的对象
saveRDS(cds, file="cds.RDS")