我正在研究一种方法来寻找受甲基化调控的基因。该方法的第一步是计算样本中每个基因的甲基化和表达数据之间的相关系数。我应用了第一个滤波器,得到了418个基因的子集,其相关系数是有意义的(小于-0.3)。
下一步是找到每个基因甲基化和表达的b样条三次回归系数。我的数据包含在一个包含两个矩阵的列表中,它看起来像这样:
> met.spl[1:5,1:5]
AAAS ABCA7 ABCC12 ABCF1 ACN9
paciente6 0.15013343 0.01489557 0.7987748 0.02826255 0.02169665
paciente7 0.10827520 0.01215497 0.8515188 0.03378193 0.02452192
paciente8 0.12423923 0.01682172 0.4182180 0.03288906 0.02046130
paciente9 0.11779075 0.02198105 0.6101996 0.06389504 0.04574667
paciente10 0.09234602 0.01526621 0.8366319 0.02868425 0.02095470
> exp.spl[1:5,1:5]
AAAS ABCA7 ABCC12 ABCF1 ACN9
paciente1 -0.82350 -1.20725 0.6000000 -0.6783 0.64500
paciente2 -1.14075 -0.59575 0.2173333 -0.2644 0.54100
paciente3 -1.43000 -1.72750 1.0015000 -1.1413 0.98625
paciente4 -1.16650 -0.76250 0.4378333 -0.6804 -0.58650
paciente5 -0.51125 -1.10325 -0.1231667 -0.1521 0.02750
为了计算第一个基因的模型系数,我做了以下步骤:
> lm(exp.spl[,1] ~ bs(met.spl[,1]), data = Lshape.spline)$coef
(Intercept) bs(met.spl[, 1])1 bs(met.spl[, 1])2 bs(met.spl[, 1])3
-1.00616163 -0.44292576 0.08767607 -0.61237162
我的目标是生成一个新对象,该对象包含418个选定基因(在行中)的系数(4列)。
似乎一个简单的sapply
可以在这里工作
library(splines)
sapply(1:ncol(exp.spl), function(i) {
lm(exp.spl[,i] ~ bs(met.spl[,i]))$coef
})
测试
met.spl<-read.table(text="AAAS ABCA7 ABCC12 ABCF1 ACN9
paciente6 0.15013343 0.01489557 0.7987748 0.02826255 0.02169665
paciente7 0.10827520 0.01215497 0.8515188 0.03378193 0.02452192
paciente8 0.12423923 0.01682172 0.4182180 0.03288906 0.02046130
paciente9 0.11779075 0.02198105 0.6101996 0.06389504 0.04574667
paciente10 0.09234602 0.01526621 0.8366319 0.02868425 0.02095470", header=T)
exp.spl<-read.table(text="AAAS ABCA7 ABCC12 ABCF1 ACN9
paciente1 -0.82350 -1.20725 0.6000000 -0.6783 0.64500
paciente2 -1.14075 -0.59575 0.2173333 -0.2644 0.54100
paciente3 -1.43000 -1.72750 1.0015000 -1.1413 0.98625
paciente4 -1.16650 -0.76250 0.4378333 -0.6804 -0.58650
paciente5 -0.51125 -1.10325 -0.1231667 -0.1521 0.02750", header=T)