Spline regression coefficients for every pair elements in two matrices R

91 views Asked by At

I'm working in a method to find genes regulated by methylation. This method consists, in the first step, in calculate the correlation coefficient between methylation and expression data for every gene in the sample. I applied this first filter and I obtain a subset of 418 genes whose correlation coefficient is significative (smaller than -0,3).

The next step is to find the coefficients of a b-spline cubic regression for methylation and expression in every gene. My data consists in a list that contains two matrices, and it looks like this:

> 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

For to calculate the coefficients of the model for the first gene, I had done this steps:

> 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 

My objective is to generate a new object containing the coefficients (4 columns) for the 418 selected genes (in the rows).

1

There are 1 answers

0
MrFlick On BEST ANSWER

Seems like a simple sapply would work here

library(splines)
sapply(1:ncol(exp.spl), function(i) {
    lm(exp.spl[,i] ~ bs(met.spl[,i]))$coef
})

tested with

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)