多时间点连续的药物处理看转录变化经典例子
最近看到了一个多组学研究，它的Bulk RNA-Seq数据集超级多，比如: GSE148692, GSE167284, GSE180023, GSE167281, GSE167282, GSE167283, GSE167285, GSE16728 and GSE180013). ，但是这些数据集目前都是currently private and is scheduled to be released on Apr 15, 2023.
虽然如此，它的实验设计还是可圈可点，值得推荐，比如 ：Figure 3. Rapid conversion of colonocytes to enterocytes after SATB2 loss
实验设计就是 (A) A time course study of colonic mucosa after a single dose of TAM treatment in Satb2cKO mice showed rapid activation of FABP6 and downregulation of CA1 at day 2 and complete replacement of CA1+ cells by FABP6+ cells by day 6. OLFM4 and LYZ1 were not robustly activated until day 30. 如下所示 ：
实验设计
因为文章的转录组数据目前不公开，所以我们先看看美图养养眼，如下所示：
主成分分析和热图
在正文重点展示了持续上升和持续下降的趋势基因集，并且进行了生物学功能数据库注释。
当然了，其它变化趋势（比如先上升再下降或）的基因集也可以进行生物学功能数据库注释，在文章的附件里面展示：
I. RNA-seq profiles and GeneOntology analysis of clusters 6-9 in SATB2 timed deletion study.
这个文章发表于2021年9月27日，美国康奈尔医学院周乔课题组在***Cell Stem Cell*** 期刊，文章标题是：《SATB2 preserves colon stem cell identity and mediates ileum-colon conversion via enhancer remodeling》，在线阅读链接 是：https://doi.org/10.1016/j.stem.2021.09.004，感兴趣的可以细读。
使用R包Mfuzz进行时间序列数据分析是很简单的
通常R包Mfuzz，Mfuzz采用了一种新的聚类算法fuzzy c-means algorithm。在文献中称这种聚类算法为soft clustering算法，相比K-means等hard clustering算法，一定程度上降低了噪声对聚类结果的干扰，而且这种算法有效的定义了基因和cluster之间的关系，即基因是否属于某个cluster, 对应的值为memebership。对于分析而言，我们只需要提供基因表达量的数据就可以了，需要注意的是，Mfuzz默认你提供的数据是归一化之后的表达量，这意味着表达量必须可以直接在样本间进行比较，对于FPKM, TPM这两种定量方式而言，是可以直接在样本间进行比较的；但是对于count的定量结果，我们必须先进行归一化，可以使用edgeR或者DESeq先得到归一化之后的数据在进行后续分析。
示例代码如下所示，其中 DEGs_exp_averp 矩阵的处理有点技巧，这里略 ：
# 首先安装R包并加载
BiocManager::install("Mfuzz")
library(Mfuzz)
## 1. 预处理：去除表达量太低或者在不同时间点间变化太小的基因等步骤
# Mfuzz聚类时要求是一个ExpressionSet类型的对象，所以需要先用表达量构建这样一个对象。
eset <- new("ExpressionSet",exprs = DEGs_exp_averp)
# 根据标准差去除样本间差异太小的基因
eset <- filter.std(eset,min.std=0)
## 2. 标准化：聚类时需要用一个数值来表征不同基因间的距离，Mfuzz中采用的是欧式距离，
# 由于普通欧式距离的定义没有考虑不同维度间量纲的不同，所以需要先进行标准化
eset <- standardise(eset)
## 3. 聚类：Mfuzz中的聚类算法需要提供两个参数，第一个参数为希望最终得到的聚类的个数，这个参数由我们直接指定；
# 第二个参数称之为fuzzifier值，用小写字母m表示，可以通过函数评估一个最佳取值
c <- 6 # 聚类个数，文章中用的6个聚类，我们也用6个
m <- mestimate(eset) # 评估出最佳的m值
cl <- mfuzz(eset, c = c, m = m) # 聚类
# 在cl这个对象中就保存了聚类的完整结果，对于这个对象的常见操作如下
cl$size # 查看每个cluster中的基因个数，看的出来与文章每个类别基因个数差了一些
[1] 2269 1982 2289 1653 1393 1729
cl$cluster[cl$cluster == 1] # 提取某个cluster下的基因
cl$membership # 查看基因和cluster之间的membership
## 4. 可视化
library(RColorBrewer)
color.2 <- colorRampPalette(rev(c("#ff0000", "Yellow", "OliveDrab1")))(1000)
mfuzz.plot(eset,cl,mfrow=c(2,3),new.window= FALSE,time.labels=colnames(DEGs_exp_averp),colo = color.2)
时间序列设计也可以结合单细胞
我看到2020的文章：《Identification of EMT signaling cross-talk and gene regulatory networks by single-cell RNA sequencing》，就是一个时间序列分析，6组不同时间点的单细胞转录组数据：
6组不同时间点的单细胞转录组数据
从标题就可以看到本文更新 EMT 过程，它涉及到如下所示的3种基因；
- epithelial markers (CDH1, EpCAM, and S100A9),
- mesenchymal markers (CDH2, FN1, and FAP),
- EMT transcription factors (TGFB1, SNAI2, and S100A6)
也是重点关心，随着时间变化的关键变量或者说指标：
随着时间变化的关键变量
学徒作业
多时间点多药物多浓度处理的多种细胞系的表达量的趋势分析；
数据集是https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE116439，大家可以任选一个细胞系的一种药物的不同浓度不同时间段处理数据，比如：
GSM3232619 A549_cisplatin_0nM_24h
GSM3232620 A549_cisplatin_0nM_2h
GSM3232621 A549_cisplatin_0nM_6h
GSM3232622 A549_cisplatin_15000nM_24h
GSM3232623 A549_cisplatin_15000nM_2h
GSM3232624 A549_cisplatin_15000nM_6h
GSM3232625 A549_cisplatin_3000nM_24h
GSM3232626 A549_cisplatin_3000nM_2h
GSM3232627 A549_cisplatin_3000nM_6h
做出趋势热图，折线图，等等！