前言

继上次我们聊完maSigpro是如何挑选不同条件相比较下的差异基因以后，这次我们来聊一下maSigpro是如何根据基因随时间的表达特征来划分cluster的，阅读本推送前请先阅读：maSigpro建模代码分析之T.fit与get.siggenes

关于see.genes

以hclust聚类为例子

我们依然用测试数据为例，并查看一下see.genes的源代码：

#读入数据
data = sigs$sig.genes$ColdvsControl$sig.profiles
edesign = edesign.abiotic

time.col = 1
repl.col = 2, 
group.cols = c(3:ncol(edesign))

#判断数据类型
time = edesign[, 1]
repvect = edesign[, 2]
groups = edesign[, 3:ncol(edesign)]
narrays <- length(time)
if (!is.null(dim(data))) {
    dat <- as.data.frame(data)
    clusterdata <- data
}
else {
    clusterdata <- data[[1]]
    dat <- as.data.frame(data$sig.profiles)
}
clusterdata <- clusterdata
if (nrow(dat) > 1) {
    dat <- as.data.frame(dat[, (ncol(dat) - length(time) + 
        1):ncol(dat)])
    count.na <- function(x) length(x[is.na(x)])
    NAs <- apply(as.matrix(dat), 1, count.na)
    count.noNa <- function(x) (length(x) - length(x[is.na(x)]))
    dat <- dat[which(apply(as.matrix(dat), 1, count.noNa) >= 
            2), ]
}
else {
    NAs <- 1
}

kdata <- NULL
out <- TRUE

我们所用的data是以Cold组与Control组相比较下的差异基因为例，一共226个差异基因，接下来我们以hclust为例，来分析聚类与时间序列的关系

k = 9
cluster.method = "hclust"
distance = "cor"

if (!is.null(clusterdata)) {
            k <- min(k, nrow(dat), na.rm = TRUE)
            if (cluster.method == "hclust") {
                if (distance == "cor") {
                  dcorrel <- matrix(rep(1, nrow(clusterdata)^2), 
                    nrow(clusterdata), nrow(clusterdata)) - cor(t(clusterdata), 
                    use = "pairwise.complete.obs")
                  clust <- hclust(as.dist(dcorrel), method = 'ward.D')
                  c.algo.used = paste("hclust", "cor", 
                                      'ward.D', sep = "_")
                }
                else {
                  clust <- hclust(dist(clusterdata, method = distance), 
                    method = agglo.method)
                  c.algo.used = paste(cluster.method, distance, 
                    agglo.method, sep = "_")
                }
                cut <- cutree(clust, k = k)
            }

首先clusterdata是指Cold组与Control组相比较下的差异基因所构成的表达矩阵，并且计算基因间的相关性，构成一个相关性矩阵，最后用全为1的矩阵（该矩阵行列数目均为差异基因的数目，223）减去基因间的相关性矩阵，从而构成dcorrel矩阵：

dcorrel <- matrix(rep(1, nrow(clusterdata)^2), 
                    nrow(clusterdata), nrow(clusterdata)) - cor(t(clusterdata), 
                    use = "pairwise.complete.obs")

最后以dcorrel矩阵作为基因的距离做系统聚类（hclust），因此到这里似乎和时间序列没有关系，只与差异基因的表达模式有关，只不过下面的图的横坐标表示所有不同的组别（包括control，cold，heat，salt组的不同处理时间点），因此我们可以看出这些差异基因在不同组别中的表达情况

这幅图区分的仅只是Cold组与Control组相比较下的差异基因的表达模式，而并不涉及到时间序列

那么这幅图代表的就是基因表达水平随时间的变化

想要弄清楚这幅图的原理，我们需要了解PlotGroups的工作原理

PlotGroups的工作原理

假设我以中位数为特征值做下面的分析

data = dat[cut == 2, ]
# cut＝2 代表取hclust聚类后的第2个cluster
summary.mode = "median"
#repvect代表生物学重复那一列
repvect = edesign[, 2]

if (summary.mode == "median") {
            ##yy代表所有列sample的所有基因表达量的中位数
            yy <- apply(as.matrix(data), 2, median, na.rm = TRUE)
            sub <- paste("Median profile of ", nrow(data), " genes")
        }

ncol = ncol(groups)
codeg <- as.character(colnames(groups))

#repvect代表生物学重复的设计矩阵
reps <- i.rank(repvect)
y <- vector(mode = "numeric", length = length(unique(reps)))
x <- vector(mode = "numeric", length = length(unique(reps)))
g <- matrix(nrow = length(unique(reps)), ncol = ncol)
#yy代表所有列sample的所有基因表达量的中位数
for (k in 1:length(y)) {
     ## y代表对相同处理的不同生物学重复的中位数表达量求均值
     y[k] <- mean(yy[reps == k], na.rm = TRUE)
     x[k] <- mean(time[reps == k])
     for (j in 1:ncol) {
        g[k, j] <- mean(groups[reps == k, j])
     }

如下面示意图所示的求法（该例子中的全部基因指的是Cold组与Control组相比较下的差异基因）：

求解画图的特征值

step1求解不同颜色的点坐标，在图上表示为不同颜色的三个点
step2求解的是拟合直线在每个时间点上的点坐标，表现为step1中每一列（每一个sample）求得的表达中位数按照生物学重复取均值

reps生物学重复

最后我们把图画出来

#导入必要的参数
ylim = NULL
xlab = "Time"
ylab = "Expression value"
cexlab = 0.8
main = NULL
sub = NULL
cex.xaxis = 1
show.fit = T
dis = NULL
min.obs = 2
alfa = 0.05
nvar.correction = FALSE
step.method = "backward"
betas = NULL
show.lines = T

if (is.null(ylim)) 
        ylim = c(min(as.numeric(yy), na.rm = TRUE), max(as.numeric(yy), 
            na.rm = TRUE))
    abcissa <- x
    xlim = c(min(abcissa, na.rm = TRUE), max(abcissa, na.rm = TRUE) * 
        1.3)
    color1 <- as.numeric(sort(factor(colnames(groups)))) + 1
    color2 <- groups
    for (j in 1:ncol) {
        color2[, j] <- color2[, j] * j
    }
    color2 <- as.vector(apply(color2, 1, sum) + 1)
    plot(x = time, y = yy, pch = 21, xlab = xlab, ylab = ylab, 
        xaxt = "n", main = main, sub = sub, ylim = ylim, xlim = xlim, 
        cex = cexlab, col = color2)
    axis(1, at = unique(abcissa), labels = unique(abcissa), cex.axis = cex.xaxis)
    if (show.fit) {
        rm <- matrix(yy, nrow = 1, ncol = length(yy))
        rownames(rm) <- c("ratio medio")
        colnames(rm) <- rownames(dis)
        fit.y <- T.fit(rm,step.method, 
            min.obs, alfa, nvar.correction)
        betas <- fit.y$coefficients
    }
    for (i in 1:ncol(groups)) {
        group <- g[, i]
        if ((show.fit) && !is.null(betas)) {
            li <- c(2:6)
            a <- reg.coeffs(coefficients = betas, groups.vector =groups.vector, 
                group = colnames(groups)[i])
            a <- c(a, rep(0, (7 - length(a))))
            curve(a[1] + a[2] * x + a[3] * (x^2) + a[4] * (x^3) + 
                a[5] * (x^4) + a[6] * (x^5) + a[7] * (x^5), from = min(time), to = max(time), 
                col = color1[i], add = TRUE, lty = li[i])
        }
        if (show.lines) {
            lx <- abcissa[group != 0]
            #y代表对相同处理的不同生物学重复的中位数表达量求均值
            ly <- y[group != 0]
            ord <- order(lx)
            lxo <- lx[ord]
            lyo <- ly[ord]
            lines(lxo, lyo, col = color1[i])
        }
    }

由此可见，每个时间点不同颜色的点代表不同的处理下的所有基因（该例子中的所有基因指的是Cold组与Control组相比较下的差异基因）表达量的中位数，当然每个点都是经过差异基因筛选（这里的差异基因筛选指的是Cold组与Control组相比较下的差异基因），然后经过hclust聚类并选取其中一类来做计算的；而每条线代表的是相同处理（相同颜色）所有基因表达量的中位数按照生物学重复计算的均值并在每个时间点之间连线而来的

值得注意的是：

if (show.fit) {
        #yy代表所有列sample的所有基因表达量的中位数
        rm <- matrix(yy, nrow = 1, ncol = length(yy))
        rownames(rm) <- c("ratio medio")
        colnames(rm) <- rownames(dis)
        fit.y <- T.fit(rm, design = dis, step.method = step.method, 
                       min.obs = min.obs, alfa = alfa, 
                       nvar.correction = nvar.correction)
        betas <- fit.y$coefficients
    }
    for (i in 1:ncol(groups)) {
        group <- g[, i]
        if ((show.fit) && !is.null(betas)) {
            li <- c(2:6)
            a <- reg.coeffs(coefficients = betas, groups.vector =groups.vector, 
                group = colnames(groups)[i])
            a <- c(a, rep(0, (7 - length(a))))
            curve(a[1] + a[2] * x + a[3] * (x^2) + a[4] * (x^3) + 
                a[5] * (x^4) + a[6] * (x^5) + a[7] * (x^5), from = min(time), to = max(time), 
                col = color1[i], add = TRUE, lty = li[i])
        }

reg.coeffs做了一个所有列sample的所有基因表达量的中位数与时间的回归分析，并用curve进行拟合曲线