我想通过从数据中提取系统发育信号来转换我的基因表达数据。我使用R包MCMCglmm。
我可以将MCMCglmm应用于表达式列之一:

require(ape)
library("MCMCglmm")
expr1 <- c(5,6,5, 11,12,13, 32,33,36)
expr2 <- c(1100,1212,1333, 32,33,36, 34, 38, 49)
expr3 <- c(32,33,36, 110,120,130, 320,330,360)
animal <- seq.int(9)
popGroup <- c(rep('A', 3),rep('B', 3), rep('C', 3))
data <- as.data.frame(cbind(expr1, expr2, expr3, animal, popGroup))
class(data$expr1)<-'integer'
class(data$expr3)<-'integer'
class(data$expr2)<-'integer'

# tree file content:
# (((1:2.0,(2:1.0,3:1.0):1.0):3.0,((4:1.0,5:1.0):1.0,6:2.0):3.0):3.0,(7:2.0,(8:1.0,9:1.0):1.0):6.0);
tree <- read.tree("tree.nwk")

prior<-list(R=list(V=1, nu=1), G=list(G1=list(V=1, nu=1)))
summary(MCMCglmm(expr1~popGroup-1, random=~animal, pedigree=tree, data=data, family="poisson", prior = prior))

但是我有超过20000这样的列。因此,我的想法是遍历所有这些对象:
for (i in 1:3) {
  M <- ( (colnames(data)[i]~popGroup-1, random=~animal, pedigree=tree, data=data, family="poisson", prior = prior))
  summM <- summary(M)
  statM <- summM$statistics[,1:2]
  print(statM)
}

问题在于在循环中定义响应变量。我尝试了很多方法,但没有任何效果。

最佳答案

这是来自MCMCglmm包的作者Jarrod Hadfield的解决方案:

不幸的是,MCMCglmm在工作时与glm略有不同
定义响应变量。您可以做的是:


固定
MCMCglmm(fixed = fixed,...)


该脚本现在可以使用了。

关于r - R:如何制作MCMCglmm循环,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/24477615/

10-10 21:55