Question: Limma for cell components effect on response
0
gravatar for antgomo
3 months ago by
antgomo0
antgomo0 wrote:

I am using limma to find the cell composition effect on response against treatment

my phenodata is as follows

...

Code  RESP Age SEX         CD8T      CD4T      NK          Bcell      Mono     Neu

30105    R    39   F    0.08900827 0.10813018 0.04065231 0.03384039 0.08240881 0.6639053
30106   NR   34   F    0.07089437 0.13440019 0.05485091 0.03168903 0.09493677 0.6452768
30119    R    73   F    0.10066220 0.23214008 0.07688516 0.04285787 0.07616766 0.4949032
30121    R    58   F    0.09589028 0.12685706 0.05535219 0.03765947 0.06127022 0.6421821
30122    R    47   F    0.04024961 0.07496977 0.02534626 0.02226978 0.05101141 0.8014216
30125   NR   66   F    0.02996649 0.05638210 0.02400648 0.02612157 0.02844631 0.8519205
30126    R    53   F    0.05369147 0.16694206 0.02350887 0.04463168 0.07090899 0.6591133
30128   NR   76   F    0.05227852 0.25069129 0.03144042 0.03237537 0.13886930 0.5275622
30134   NR   47   F    0.08675013 0.17954926 0.03897045 0.08915519 0.10016838 0.5315112
30135    R    55   F    0.06359675 0.15270431 0.03647699 0.04048208 0.07990695 0.6537159

...

So i was thinking i something like that, let's say that i want to test each effect separately to get the effect of this particular cell component in response to treatment

design <- model.matrix(~0+pd$RESP+pd$RESP:pd$Neu+ pd$SEX + pd$Age)

In this case, i have to test the5 and 6 coefficients, am i right?

pd$RESPNR pd$RESPR pd$SEXM pd$Age pd$RESPNR:pd$Neu pd$RESPR:pd$Neu
  1         0        1                      0     39        0.0000000       0.6639053
  2         1        0                      0     34        0.6452768       0.0000000
  3         0        1                      0     73        0.0000000       0.4949032
  4         0        1                      0     58        0.0000000       0.6421821
  5         0        1                      0     47        0.0000000       0.8014216
  6         0        1                      0     66        0.0000000       0.8519205

However, i don't know if this is the real design that i have to follow, maybe i have to consider splines and hence

X <- ns(df$Neu, df=3)

design <- model.matrix(~pd$RESP*X+ pd$Sex + pd$Age)

Any hint would be greatly appreciated

Many thanks in advance

ADD COMMENTlink modified 3 months ago by Aaron Lun25k • written 3 months ago by antgomo0
Answer: Limma for cell components effect on response
0
gravatar for Aaron Lun
3 months ago by
Aaron Lun25k
Cambridge, United Kingdom
Aaron Lun25k wrote:

You should throw all the cell compositions into the model at once. This protects you against co-linearity in the compositions. For example, if CD8 and CD4 T cell abundances are strongly correlated, a single model that contains terms for all cell abundances will not be able to confidently associate a particular gene to either cell type. This is the correct outcome - using separate models will suggest that a gene is associated with both CD8 and CD4 T cell abundance, which could be misleading.

It also protects against variance inflation and loss of power from a mis-specified model. For example, let's say that the expression profiles of many genes are associated with B cell abundance. But if you don't include B cells in your (say) neuron-specific model, the apparent variance for those genes will be inflated. This reduces the power to detect any association with neuron abundance - but more importantly, due to empirical Bayes shrinkage, the inflated variances for these genes will also reduce power for other genes, which is not good.

So just throw in pd$RESP:pd$XXX terms for all cell types XXX. And hope that you have enough residual d.f. to fit the model and estimate the variances. I wouldn't worry about splines unless you really have loads and loads of samples to be able to afford them.

ADD COMMENTlink modified 3 months ago • written 3 months ago by Aaron Lun25k

Many thanks Aaron,

so i will use

design <- model.matrix(~pd$RESP*cell_comp+ pd$Sex + pd$Age)

without splines, in this case, i will have to look individually for each pd$RESP:pd$XXX individually, isn't it?

if i want to check the effect of cell_comp in DE genes between R and NR, i have to do contrasts

con <- makeContrasts( RESP_XXX - NRESP_XXX,
levels=design) # replacing ":" with "_" in colnames(design)
ADD REPLYlink modified 3 months ago • written 3 months ago by antgomo0

Yes and yes. You may also want to log-transform the components, I would expect log-expression to increase linearly with log-abundance rather than linearly with raw abundance.

ADD REPLYlink written 12 weeks ago by Aaron Lun25k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 286 users visited in the last hour