2
$\begingroup$

I have FPKM values for Normal(80 cases), Group1(135 cases), Group2(147 cases), and Group3(102 cases) for 160 genes (Groups were obtained by clustering tumor data).

I want to perform t-test (two-sided, unpaired, unequal variance) between Normal and Group1, Normal and Group2, and Normal and Group3 such that I get p.value for each gene.
The t.test function in R returns single p.value and not gene-wise (I don't want to do t-test in Excel).

Which package/function should I use to do this in R?

$\endgroup$

1 Answer 1

4
$\begingroup$

You can simply apply() t.test() to your matrix. In general, though, I expect lmFit() from the limma package will be simpler to use.

Here is an example with t.test():

m = matrix(runif(100, max=c(rep(1, 50), rep(1.5, 50))), nrow=10) # make some data
apply(m, 1, function(x) t.test(x[1:5], x[6:10])$p.value)

Neither t.test() nor lmFit() should be expected to produce the same results as excel, since excel is inferior in every way and should never be relied upon for any statistical analysis (plus lmFit is doing something much much more complicated than a simple t-test).

As an aside, FPKMs are highly problematic for use in statistics. If possible, I strongly encourage you to get either per-gene counts or per-gene estimated counts and to use those in limma or another differential expression package instead, as the results will be much more reliable. I presume, though, that you're using TCGA data, so you may be stuck with FPKMs.

$\endgroup$
7
  • 1
    $\begingroup$ I totally agree with Devon's aside comment. FPKM values are not suitable for (sound) statistics and only informative for heatmaps etc. For proper statistics you should use count data, all statisticians agree on that. $\endgroup$
    – benn
    Commented Apr 15, 2018 at 9:20
  • $\begingroup$ @b.nota Source? $\endgroup$
    – SmallChess
    Commented Apr 15, 2018 at 13:02
  • $\begingroup$ @SmallChess, for example paper about edgeR, but also the more difficult stuff of Gordon Smyth and read the discussions on bioconductor support. $\endgroup$
    – benn
    Commented Apr 15, 2018 at 13:48
  • $\begingroup$ I tried t.test()Group1.ttest <- t.test(data.Normal, data.Group1, "two.sided", paired = FALSE, var.equal = FALSE)but the function does not give me p.values row-wise. The result that I obtained for above code> Group1.ttest$p.value [1] 5.632393e-48I tried lmFit() fac <- factor(c(rep(0,80),rep(1,135))) fit = lmFit(Normal.Group1, design=model.matrix(~ fac)) e = eBayes(fit) tab<-topTable(e, number=nrow(Normal.Group1)) but the p-values obtained are different from the ones obtained when I do the t.test on the same data in Excel. I tried log2(Normal.Group1) but still the results differ. $\endgroup$
    – user98059
    Commented Apr 15, 2018 at 13:55
  • 2
    $\begingroup$ The results should differ, excel should never be relied on for anything important. lmFit is also a more complicated computation than a simple t-test. Please see my updated reply. $\endgroup$
    – Devon Ryan
    Commented Apr 15, 2018 at 18:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.