背景
前面我们画过WGCNA得到模块以后,用模块包含的基因表达量做拟合曲线。
例如文章PMID: 26460053中的图片

使用场景
本质是用一条趋势线代表多个点或线的变化趋势,多用于时间序列数据。
场景一:WGCNA得到的模块内的基因的变化趋势
场景二:聚类分析得到的cluster内的基因的变化趋势
输入数据
之前我们介绍过wcgna输出各个gene module内的基因表达量,每个module保存在一个csv文件中。
例子里的数据是时间序列,作者提取表达量最高的30个基因,取平均值(averaged gene expression of top 30)绘制拟合曲线。
一次画一个module
这套输入数据里,每种细胞有5个sample,假装是5个时间点,4次重复,取前30个基因在4次重复里的平均值
specmod<-read.csv("pink.csv",row.names = 1)[1:30,]
t<-data.frame(tp=1:5,tpm_mean = log2(c(mean(as.matrix(specmod[,1:4])),mean(as.matrix(specmod[,5:8])),mean(as.matrix(specmod[,9:12])),mean(as.matrix(specmod[,13:16])),mean(as.matrix(specmod[,17:20])))))
head(t)

画图
library(ggplot2)
ggplot(t,aes(x = tp,y = tpm_mean)) +
geom_smooth(color="black") + #拟合曲线的颜色
labs(x = "Days after injury", y = "Expression Level", title = "black") +
theme_bw() + #去除背景色
theme(panel.grid =element_blank()) + #去除网格线
theme(panel.border = element_blank()) + #去除外层边框
theme(axis.line = element_line(colour = "grey")) + #坐标轴画成灰色
theme(axis.ticks = element_blank()) + #取掉坐标轴上的刻度线
geom_hline(yintercept = t$tpm_mean[1],linetype="dashed") #在第一个时间点处画虚线

默认用loess模型拟合,你还可以添加method参数,使用其他方法拟合,例如:method = "lm"
, 或 "glm", "gam", "loess", "rlm"
一次画多个module
读入当前文件夹中的多个module文件,合并。
fnames<-Sys.glob("*.csv")
fnames
fdataset<-lapply(fnames,read.csv)
names(fdataset) <- fnames
library(plyr)
result <- ldply(fdataset, data.frame)
head(result)
df<-result[,2:22]
df$module<-unlist(strsplit(result$.id,split = ".csv"))
row.names(df)<-df$X
df1<-df[2:22]
df1前几行如下:

画图
library(ggplot2)
plist=list()
tag = 0
for (mod in unique(df1$module)){
tag = tag + 1
specmod <- df1[df1$module == mod,]
t<-data.frame(tp=1:5,tpm_mean = log2( c(mean(as.matrix(specmod[1:30,1:4])),mean(as.matrix(specmod[1:30,5:8])),mean(as.matrix(specmod[1:30,9:12])),mean(as.matrix(specmod[1:30,13:16])),mean(as.matrix(specmod[1:30,17:20])))))
plist[[tag]]<-ggplot(t,aes(x = tp, y = tpm_mean)) +
geom_smooth(color = mod) + #拟合曲线的颜色
labs(x = "Days after injury", y = "Expression Level", title = mod) +
theme_bw() + #去除背景色
theme(panel.grid =element_blank()) + #去除网格线
theme(panel.border = element_blank()) + #去除外层边框
theme(axis.line = element_line(colour = "grey")) + #坐标轴画成灰色
theme(axis.ticks = element_blank()) + #取掉坐标轴上的刻度线
geom_hline(yintercept = t$tpm_mean[1],linetype="dashed") #在第一个时间点处画虚线
}
require(cowplot)
plot_grid(plotlist=plist, ncol=2)
ggsave(file="myfitting.pdf")

8个模块,各自的top30基因在时间序列上的表达情况。