普通的Bulk转录组也能做“拟时序”分析吗？

[ixxmu]

https://mp.weixin.qq.com/s/Y203ep6cpR6DyTjoMKgflA

[ixxmu]

普通的Bulk转录组也能做“拟时序”分析吗？ by BioInfo

写在前面：Debug不易，码字也不易，如果你采用了本文的代码，请引用一下下面的几篇文章，拜托了。

那么，正文开始

一般来说，拟时序分析（Pseudotime analysis）通常出现在单细胞转录组数据中，通过模拟细胞在发育或者分化过程中的时间顺序来揭示不同细胞状态之间的动态变化。在单细胞转录组数据中常用Monocle或Slingshot来进行拟时序分析。

但是普通的Bulk转录组是否也能进行拟时序分析呢？

例如在这篇两文章中，作者都是利用Bulk转录组数据做了拟时序分析，对基因表达的先后进行排序，从而使得基因拥有了“时序性”。

该文章阐述了六倍体小麦胚胎发育过程中的动态染色质调控程序。在Fig5中，作者利用Bulk的RNA-seq和ATAC-seq数据进行了“拟时序分析”为基因赋予了“时序性”，从而搭建了具有“时序性”的转录调控网络。

同样地，在第二篇文章中，作者则通过类似的方式搭建了小麦再生过程中的“时序性”的转录调控网络。

两篇文章中使用的方法均自之前的一篇文章，利用拟时序分析来对玉米和大麦花序发育过程中的转录组进行比较。

但其实，这种利用bulk转录组做的拟时序分析与单细胞转录组中的拟时序分析是完全不同的。

单细胞转录组中的拟时序分析是对细胞进行分析，根据每个细胞的基因表达谱数据，为每个细胞分配发育时间点的信息。例如一个样本中同时存在祖先细胞和分化后的细胞，如果将发育最早的祖先细胞的时间设置为0，将发育最晚的分化后的细胞的时间设置为10，那么拟时序分析的作用就是为这个样本中所有细胞都分配一个在0-10之间的发育时间点，越接近0就说明这个细胞与祖先细胞越相似，而越接近于10则说明这个细胞的分化程度越高。

但是在bulk转录组中，拟时序的对象不再是细胞，而是基因（严格意义上来说，不清楚这种分析是否也能被称为拟时序分析，在本文中暂且将其称为拟时序分析）。在bulk转录组中，一般只有少数样本，最多可能20-30个样本点，其分辨率远远达不到单细胞转录组中的几千甚至上万个细胞的分辨率（在单细胞转录组中，1个细胞即可被认为是一个样本点）。因此用bulk转录组进行拟时序分析时，首先需要将少数几个样本点变成500或者更多个样本点，这样才能更好地根据基因表达的峰值的时间点来对基因表达的先后顺序进行排序，进而得到“拟时序”的基因。

具体做法如下：

假设我们有5个样本点（A, B, C, D, E）的1000个基因表达量的数据，我们希望对这1000个基因进行所谓的“拟时序”分析，以获得这1000个基因表达的先后顺序。

Step1：对5个样本点进行PCA分析，并通过PC1和PC2计算样本点之间的欧式距离获得时间线，最终将时间线缩放到0-10之间。

library(splines)
library(factoextra)
library(FactoMineR)
library(scales)
#读取基因表达的TPM矩阵
tpm=read.table("gene.TPM.txt",header=T,row.names=1)
#对5个样本进行PCA分析
pca=PCA(t(tpm),scale.unit=T,ncp=5,graph=F)
#提取PC1和PC2的结果，计算5个样本点的欧式距离
res=pca$ind$coord
res=res[,1:2]
dis=dist(res)
dis=as.matrix(dis)
#利用相邻样本点的欧式距离计算5个样本点原始的时间线
point1=0
point2=point1+dis[1,2]
point3=point2+dis[2,3]
point4=point3+dis[3,4]
point5=point4+dis[4,5]
raw.timeline=c(point1,point2,point3,point4,point5)
#通过scales包中的rescale函数将原始的时间线经过线性变换为0-10之间的新的时间线
new.timeline=rescale(raw.timeline,to=c(0,10))

Step2：利用modelr中的loess函数，根据新时间线，产生500个时间点的表达值。

library(tidyverse)
library(modelr)
#对tpm值进行zscale标准化
tpm.z=as.data.frame(t(apply(tpm,1,scale)))
names(tpm.z)=names(tpm)
#将5个时间点的zcale之后的表达值变成500个点的表达值
grid=data.frame(time=seq(0,10,length.out=500))
pseudotime.model.fun=function(value){time=new.timeline;data = tibble(value=value,time =time);model=loess(value ~ time,data);predict=grid %>% add_predictions(model);return(predict)}
res=apply(tpm.z,1,pseudotime.model.fun)
results=res %>% reduce(bind_cols)
exp500=as.data.frame(t(results[,seq(0,ncol(results),2)]))
names(exp500)=results$time...1
row.names(exp500)=row.names(tpm)

Step3: 利用atan2函数对基因表达进行排序，定义基因表达的先后顺序

library(ComplexHeatmap)
library(factoextra)
library(FactoMineR)
#对基因进行PCA分析，得到每个基因的Dim1和Dim2
gene.pca=PCA(exp500,scale.unit=T,ncp=5,graph=F)
res=gene.pca$ind$coord
res=res[,1:2]
exp500=cbind(exp500,res)
#用基因的Dim1和Dim2，进行排列组合，计算四种atan2值
exp500$atan2.1=atan2(exp500$Dim.1,exp500$Dim.2)
exp500$atan2.2=atan2(exp500$Dim.1,-exp500$Dim.2)
exp500$atan2.3=atan2(-exp500$Dim.1,exp500$Dim.2)
exp500$atan2.4=atan2(-exp500$Dim.1,-exp500$Dim.2)
#分别根据四种atan2值从低到高进行排序，然后用ComplexHeatmap可视化结果选定排序方式
order1=arrange(exp500,exp500$atan2.1)
order2=arrange(exp500,exp500$atan2.2)
order3=arrange(exp500,exp500$atan2.3)
order4=arrange(exp500,exp500$atan2.4)
p1=Heatmap(order1[,1:500],cluster_rows = F,cluster_columns = F,show_row_names = F,show_column_names = F,column_title = "Order1",heatmap_legend_param = list(title="Order1",legend_height=unit(2,"cm")),col=colorRampPalette(rev(brewer.pal(n = 11, name="RdYlBu")))(100))
p2=Heatmap(order2[,1:500],cluster_rows = F,cluster_columns = F,show_row_names = F,show_column_names = F,column_title = "Order2",heatmap_legend_param = list(title="Order2",legend_height=unit(2,"cm")),col=colorRampPalette(rev(brewer.pal(n = 11, name="RdYlBu")))(100))
p3=Heatmap(order3[,1:500],cluster_rows = F,cluster_columns = F,show_row_names = F,show_column_names = F,column_title = "Order3",heatmap_legend_param = list(title="Order3",legend_height=unit(2,"cm")),col=colorRampPalette(rev(brewer.pal(n = 11, name="RdYlBu")))(100))
p4=Heatmap(order4[,1:500],cluster_rows = F,cluster_columns = F,show_row_names = F,show_column_names = F,column_title = "Order4",heatmap_legend_param = list(title="Order4",legend_height=unit(2,"cm")),col=colorRampPalette(rev(brewer.pal(n = 11, name="RdYlBu")))(100))
p1+p2+p3+p4

从结果中可以看到，Order3是符合预期的，因此保留Order3中的结果。

PS：Debug不易，码字也不易，如果你采用了本文的代码，请引用一下上述的几篇文章。

Zhao L, Yang Y, Chen J, et al. Dynamic chromatin regulatory programs during embryogenesis of hexaploid wheat. Genome Biol. 2023;24(1):7. Published 2023 Jan 13. doi:10.1186/s13059-022-02844-2

Liu X, Bie XM, Lin X, et al. Uncovering the transcriptional regulatory network involved in boosting wheat regeneration and transformation [published correction appears in Nat Plants. 2024 Jan;10(1):194. doi: 10.1038/s41477-024-01619-w]. Nat Plants. 2023;9(6):908-925. doi:10.1038/s41477-023-01406-z

Leiboff S, Hake S. Reconstructing the Transcriptional Ontogeny of Maize and Sorghum Supports an Inverse Hourglass Model of Inflorescence Development. Curr Biol. 2019;29(20):3410-3419.e3. doi:10.1016/j.cub.2019.08.044