前言

上一次我们说到了make.design.matrix，maSigpro建模代码分析之make.design.matrix
是如何构建建模所用的设计矩阵，那么这次我们将讨论下如何计算模型的p值，首先是计算模型回归系数总体显著性的函数p.vector()

基于上一个函数make.design.matrix我们得到了设计矩阵：

design <- make.design.matrix(edesign, degree = 3)

其中：

我们的对象design里面包含dis矩阵

并且我们的基因表达矩阵data为：

data

p.vector

首先，p.vector函数包括如下几个参数

p.vector <- function (data, design = NULL, Q = 0.05, MT.adjust = "BH", min.obs = 3, counts=FALSE, family=NULL, theta=10, epsilon=0.00001) 

显然我们的对象design是一个list：

design

在接下来的例子中，我们以”高斯分布“举例：

if (is.list(design)) {
  dis <- as.data.frame(design$dis)
  groups.vector <- design$groups.vector
  edesign <- design$edesign
}

#判断是否为高斯分布
if (is.null(family)) {
  if(!counts) { family=gaussian() }
  if(counts)  { family=negative.binomial(theta) }
}

#获得基因表达矩阵数据
min.obs = 3
dat <- as.matrix(data)
#去除NA值
dat <- na.omit(dat)
dat <- dat[, as.character(rownames(dis))]
G <- nrow(dat)
#这一步感觉像过滤表达量低的基因
count.na <- function(x) (length(x) - length(x[is.na(x)]))
dat <- dat[apply(dat, 1, count.na) >= min.obs, ]
#一共有1000个基因，g = 1000
g <- dim(dat)[1]
n <- dim(dat)[2]
p <- dim(dis)[2]
p.vector <- vector(mode = "numeric", length = g)

获得必要的基本信息以后就可以开始计算p值了，那么我们以高斯分布为例：

##对每一个基因进行多项式回归建模
for (i in 1:g) {
  # y为各个基因表达量
  y <- as.numeric(dat[i, ])
  div <- c(1:round(g/100)) * 100
  #由设计矩阵参数所得的多项式回归模型
  model.glm<- glm(y~.,data=dis , family="gaussian", epsilon=0.00001)
  if(model.glm$null.deviance==0){ 
    p.vector[i]=1 
  } else{
    # "1"(model.glm.0)代表的是对照组模型
    model.glm.0<-glm(y~1, family="gaussian", epsilon=0.00001)
    #计算p值
    test<-anova(model.glm.0,model.glm,test="F")
      if (is.na(test[6][2,1])) {
        p.vector[i]=1
      } else{
        p.vector[i]=test[6][2,1]
      }
   }
}

model.glm
其中model.glm<- glm(y~.,data=dis , family="gaussian",epsilon=0.00001)代表的是多项式回归方程模型，里面有各个多项式的回归系数 model.glm.0
其中model.glm.0<-glm(y~1, family="gaussian", epsilon=0.00001)代表的是原假设模型，即各个回归系数均为0 test

这里的p.vector代表的是两个模型model.glm和model.glm.0利用anova函数做F检验的结果，即检验多项式回归方程的各个回归系数总体显著性
之后根据p值来做矫正：

p.adjusted <- p.adjust(p.vector, method = "BH", n = length(p.vector))
#筛选出总体回归系数显著的基因
genes.selected <- rownames(dat)[which(p.adjusted <= 0.05)]

FDR <- sort(p.vector)[length(genes.selected)]

#SELEC即为我们筛选出的总体回归系数显著的基因
SELEC <- as.matrix(as.data.frame(dat)[genes.selected, ])
if (nrow(SELEC) == 0) 
  print("no significant genes")

p.vector <- as.matrix(p.vector)
rownames(p.vector) <- rownames(dat)
colnames(p.vector) <- c("p.value")

output <- list(SELEC, p.vector, p.adjusted, G, g, FDR, 
               nrow(SELEC), dis, dat, min.obs, 0.05, groups.vector, edesign, family)
names(output) <- c("SELEC", "p.vector", "p.adjusted", "G", 
                   "g", "FDR", "i", "dis", "dat", "min.obs", "Q", "groups.vector", 
                   "edesign", "family")
output

SELEC即为我们筛选出的总体回归系数显著的基因
那么最终的output主要是关于p值的相关的一些数据，这就是p.vector函数的主要作用，后续会根据p值对基因做一些筛选

下期预告：maSigpro建模代码分析之T.fit