Difficulty Reproducing log Fold Change Calculation Using tmod Package in R

33 views Asked by At

I want to calculate the minimum significant differences (MSD) of genes expressed in treated cells compared to control cells. The MSD can be used to rank genes for gene set enrichment analysis. Therefore, I have a dataframe containing protein abundance values of different genes per row. Columns are the various replicates and treatments I found the package tmod, that can perform the calculation of the MSD. Tmod is based on the limma package. Besides the MSD it performes various other calculations, i.e. the log fold change. However, I am not able to reproduce the calculation of the logFC as performed by tmod.

Here is an example to reproduce the problem:

> library(limma)
> library(tmod)
> numerator_columns <- c("ERU_4", "ERU_5", "ERU_6")
> denominator_columns <- c("C_1", "C_2", "C_3")
> test_df <- data.frame(
+   "C_1" = c(8, 7, 3),
+   "C_2" = c(6, 8, 9),
+   "C_3" = c(1, 2, 7),
+   "ERU_4" = c(3, 1, 10),
+   "ERU_5" = c(2, 2, 1),
+   "ERU_6" = c(10, 10, 8)
+ )
> all_cols <- colnames(test_df)
> cond <- rep(-1, length(all_cols))
> cond[all_cols %in% numerator_columns] <- "Num"  
> cond[all_cols %in% denominator_columns] <- "Denom"  
> cond <- as.factor(cond)
> contr <- c(
+   "Num-Denom"
+ )
> design <- model.matrix(~0+cond)
> colnames(design) <- levels(cond)
> contrast =  makeContrasts(contrasts= contr,levels=design)
> fit1 <- lmFit(test_df, design)
> fit2 <- contrasts.fit(fit1,contrasts = contrast)
> fit3 <- eBayes(fit2)
> result <- tmodLimmaTopTable(fit3, adjust.method = "BH", confint = 0.95)
> print(result)

Tmod produces the following dataframe:

> print(result)
  logFC.Num-Denom   t.Num-Denom msd.Num-Denom SE.Num-Denom   d.Num-Denom ciL.Num-Denom ciR.Num-Denom qval.Num-Denom
1    8.881784e-16  2.687021e-16     -7.201933     3.305439  3.595630e-15     -7.201933      7.201933              1
2   -1.333333e+00 -4.033756e-01     -5.868599     3.305439 -5.397759e+00     -8.535266      5.868599              1
3    0.000000e+00  0.000000e+00     -7.201933     3.305439  0.000000e+00     -7.201933      7.201933              1

However, I would expect the logFC to look like:

> print(result)
  logFC.Num-Denom   t.Num-Denom msd.Num-Denom SE.Num-Denom   d.Num-Denom ciL.Num-Denom ciR.Num-Denom qval.Num-Denom
1    0.000000e+00  2.687021e-16     -7.201933     3.305439  3.595630e-15     -7.201933      7.201933              1
2   -3.870231e-01 -4.033756e-01     -5.868599     3.305439 -5.397759e+00     -8.535266      5.868599              1
3    0.000000e+00  0.000000e+00     -7.201933     3.305439  0.000000e+00     -7.201933      7.201933              1

I also tried to change the way contrasts are formed:

> cont <- c(
+   "Num/Denom"
+ )

# or

> cont <- c(
+   "10^(log10(Num)-log10(Denom))"
+ )

However, this only produces NA values or Inf.

Can you help me to understand how tmod calculates the logFC and the MSD? Do you know different ways to calculate the MSD?

0

There are 0 answers