I ran a DGE analysis with 4 different contrasts with lmFit() + contrasts.fit() in limma, and I'm interested in the overlap of the four different contrasts. I identified the genes that are indeed differentially expressed in all four contrasts, and coded that in tfit$genes$test ("yes" if differentially expressed, "no" if not). The MArrayLM object is pasted below.
I did manage to run a Kegg pathway enrichment analysis for each of the contrasts individually, but I'm after a way to run a single analysis for the set of DE genes in all contrasts. Is that feasible?
An object of class "MArrayLM" $coefficients Contrasts ARvsCA ARvsEU ARvsAU ARvsNZ 1 -0.18067896 -0.2044603 -0.22881771 -0.1833862 2 -0.04079345 -1.1859285 -0.39206763 -0.3653143 3 -0.12733594 0.1763288 -0.07934863 -0.1252855 4 0.07648264 0.6875827 0.13266024 0.3442508 5 0.09678434 0.4514540 0.13137207 0.2943875 10387 more rows ... $stdev.unscaled Contrasts ARvsCA ARvsEU ARvsAU ARvsNZ 1 0.1820675 0.2300422 0.1936519 0.1780754 2 0.1499697 0.1868217 0.1574127 0.1458230 3 0.1480576 0.1898012 0.1567854 0.1442379 4 0.1883755 0.2695961 0.2030905 0.1910940 5 0.1406236 0.1774499 0.1479352 0.1374007 10387 more rows ... $sigma [1] 1.2656682 1.3338325 0.5922579 1.9853202 0.9379173 10387 more elements ... $df.residual [1] 21 21 21 21 21 10387 more elements ... $cov.coefficients Contrasts Contrasts ARvsCA ARvsEU ARvsAU ARvsNZ ARvsCA 0.3666667 0.2000000 0.2 0.2000000 ARvsEU 0.2000000 0.5333333 0.2 0.2000000 ARvsAU 0.2000000 0.2000000 0.4 0.2000000 ARvsNZ 0.2000000 0.2000000 0.2 0.3428571 $rank [1] 5 $genes gene_id line test 1 gene12245 1 no 2 gene12244 2 no 3 gene12247 3 no 4 gene12246 4 no 5 gene12241 5 no 10387 more rows ... $Amean 1 2 3 4 5 4.749884 8.665232 5.995541 4.329715 6.534481 10387 more elements ... $method [1] "ls" $design AR_heads AU_heads CA_heads EU_heads NZ_heads 1 1 0 0 0 0 2 1 0 0 0 0 3 1 0 0 0 0 4 1 0 0 0 0 5 1 0 0 0 0 21 more rows ... $contrasts Contrasts Levels ARvsCA ARvsEU ARvsAU ARvsNZ AR_heads 1 1 1 1 AU_heads 0 0 -1 0 CA_heads -1 0 0 0 EU_heads 0 -1 0 0 NZ_heads 0 0 0 -1 $df.prior [1] 2.659777 $s2.prior [1] 0.6719199 $s2.post [1] 1.4973681 1.6546414 0.3868724 3.5739382 0.8563319 10387 more elements ... $df.total [1] 23.65978 23.65978 23.65978 23.65978 23.65978 10387 more elements ... $t Contrasts ARvsCA ARvsEU ARvsAU ARvsNZ 1 -0.1937939 -0.2378606 -0.3853473 -0.2105622 2 0.0000000 -4.3627269 -1.2572034 -1.2144966 3 0.0000000 0.3288759 0.0000000 0.0000000 4 0.0000000 1.0792898 0.0000000 0.5722939 5 0.0000000 1.9118965 0.0000000 1.2338677 10387 more rows ... $p.value Contrasts ARvsCA ARvsEU ARvsAU ARvsNZ 1 0.5071515 0.5252303689 0.4194142 0.4945336 2 0.8719553 0.0001139932 0.1181080 0.1248538 3 0.5476869 0.3794963643 0.7396859 0.5572903 4 0.8441133 0.2050344009 0.7492352 0.3837871 5 0.6637349 0.0347921880 0.5483615 0.1158894 10387 more rows ... $treat.lfc [1] 0.1375035
Yes, it's easy enough. But you need to have generally recognized gene IDs (usually Entrez Gene Ids) order to run kegga(). Your gene_ids don't seem to be Entrez Ids. Have you just anonymized them?
Hi Gordon,
Sorry I reported the wrong MArrayLM fit. Before using it as input for kegga(), I did change the gene_id into RefSeq identifiers as below:
$genes
gene_id line test
1 105678280 1 no
2 105678292 2 no
3 105678279 3 no
4 105678278 4 no
5 105678296 5 no
10387 more rows ...