Differential Gene Expressionand Pathway Enrichment: tGFR vs. cGFR

Pipeline overview

This script performs pairwise testing for gene-expression results via Affymetrix gene chip, visualizes the results, places the differentially-expressed genes into pathways via GO and Reactome, and performs pathway enrichment analyses.

Sources: The second half of this script (pathway enrichment ) draws from https://bioconductor.org/packages/devel/workflows/vignettes/maEndToEnd/inst/doc/MA-Workflow.html

Preparation for adaptation to a different pairwise comparison within the project

You should already have prepared test matrix files according to limma protocols

  1. The strings “tGFR.cGFR” and “tGFR vs. cGFR” are used throughout the script to denote the test being performed. These should be bulk find-and-replaced throughout the script
  2. The “targets” matrix in section “3. Set up test matrices” must be changed to the appropriate test
  3. Sections “5. Load in the metadata” and “6. Set up the test” must be customized to the desired test

Results files;

1. Prepare the base environment

2. Prepare Annotation Table to yield useful mappings of probe to gene

Annotation <- read.csv("HTA-2_0.na35.2.hg19.transcript.csv/HTA-2_0.na35.2.hg19.transcript.csv", header=TRUE, skip=19)
probe.gene.map <- read.csv("probe.gene.map.csv")

3. Set up test matrices

Differential abundance can only be performed on pairs of sample sets. Each pairwise test has to be configured as a separate file. Importing the test matrices:

targets <- read.csv("targets.GFRonly.txt", as.is=T)

all of the samples and metadata are listed here:

targets 
##                                      FileName Strain Treatment Experiment
## 1   20160505_HTA2.0_0318_WTH_3R_(HTA-2_0).CEL    GFR   control          1
## 2   20160505_HTA2.0_0319_WTH_4R_(HTA-2_0).CEL    GFR      drug          1
## 3 20160616_HTA2.0_0335_WTH_3_R2_(HTA-2_0).CEL    GFR   control          2
## 4 20160616_HTA2.0_0336_WTH_4_R1_(HTA-2_0).CEL    GFR      drug          2

4. Define the subset of .CEL files, read them in, and normalize them

dat <- read.celfiles(targets$FileName)
## Reading in : 20160505_HTA2.0_0318_WTH_3R_(HTA-2_0).CEL
## Reading in : 20160505_HTA2.0_0319_WTH_4R_(HTA-2_0).CEL
## Reading in : 20160616_HTA2.0_0335_WTH_3_R2_(HTA-2_0).CEL
## Reading in : 20160616_HTA2.0_0336_WTH_4_R1_(HTA-2_0).CEL
project.bgcorrect.norm.avg <- rma(dat, background=TRUE, normalize=TRUE)
## Background correcting
## Normalizing
## Calculating Expression
project.bgcorrect.norm.avg.filt <- exprs(project.bgcorrect.norm.avg)[which(rownames(project.bgcorrect.norm.avg) %in% Annotation[,1]),]

5. Load in the metadata

Exp <- factor(targets$Experiment)
Strain <- factor(targets$Strain, levels=c("GF","GFR"))
Treatment <- factor(targets$Treatment, levels=c("control","drug"))

6. Define the test

This is important decision; these chips were ran from different experimental batches, thus I have controlled for experiment, essentially running a paired analysis. In an earlier iteration, I did not control for experiment and found very few differences due to the high variance between experimental replicates

design <- model.matrix(~Exp+Treatment)

7. Run “fit” algorithms and separate out the significantly different probes

eset <- project.bgcorrect.norm.avg.filt
fit <- lmFit(eset, design)
fit <- eBayes(fit, trend=T, robust=T)
## Warning: 65 very small variances detected, have been offset away from zero
results <- decideTests(fit, p.value = 0.05, method = "separate", adjust.method = "fdr")
write.table(summary(results), "pairwise_test_results/tGFR.cGFR/summary.tGFR.cGFR.csv", quote=F)

This is a simple summary of the differences, only the far right column

summary(results)
##        (Intercept)  Exp2 Treatmentdrug
## Down             0    13            22
## NotSig          87 70484         70417
## Up           70436    26            84

8. Generate tables of the differentially expressed genes and attach gene name/function information

This script also generates a temporary differential expression table for ALL probes which can be used later for comparisons with other pairwise-test results

library(data.table)
library(dplyr)
#full table
fit.tGFR.cGFR <- topTable(fit, n = nrow(fit))
#diff table only, must be customized based on summary results
diff <- topTable(fit, p.value = 0.05,  n = Inf, coef="Treatmentdrug")
setDT(diff, keep.rownames = TRUE)[]
colnames(diff)[1] <- "probe"
diff <- dplyr::left_join(diff,probe.gene.map,by="probe")
## Warning: Column `probe` joining character vector and factor, coercing into
## character vector
write.table(diff, "pairwise_test_results/tGFR.cGFR/diff.tGFR.cGFR", quote = F)
write.table(diff$name, "pairwise_test_results/tGFR.cGFR/diff.tGFR.cGFR.gene.list", quote = F)
kable(diff) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")
probe logFC AveExpr t P.Value adj.P.Val B name description
TC10002567.hg.1 1.9922305 4.099263 15.511819 1.00e-07 0.0059623 6.331456 — —
TC11003279.hg.1 1.9618456 4.579688 15.042410 2.00e-07 0.0059623 6.200876 SYTL2 synaptotagmin-like 2
TC01005754.hg.1 2.2354093 6.161493 13.887169 4.00e-07 0.0059623 5.843220 — —
47422416_st -1.5784456 2.295995 -13.774806 4.00e-07 0.0059623 5.805414 — —
TC10002568.hg.1 2.2178789 5.768782 13.727863 5.00e-07 0.0059623 5.647471 — —
TC04002616.hg.1 1.7139339 3.857110 13.231673 5.00e-07 0.0059623 5.614363 — —
TC03000565.hg.1 2.0089888 5.876989 12.777460 7.00e-07 0.0059623 5.443360 — —
TC18000797.hg.1 1.6680374 5.184303 12.558056 8.00e-07 0.0059623 5.356845 — —
TC06000399.hg.1 1.6760826 5.256095 12.446007 9.00e-07 0.0059623 5.311635 HLA-DRA major histocompatibility complex, class II, DR alpha
TC6_dbb_hap3000083.hg.1 1.6723321 5.491333 12.185183 1.00e-06 0.0059623 5.203619 HLA-DRA major histocompatibility complex, class II, DR alpha
47422463_st -1.2844631 1.906223 -12.040285 1.10e-06 0.0059623 5.141884 — —
TC6_cox_hap2000090.hg.1 1.5840786 5.272219 11.884194 1.30e-06 0.0059623 5.073956 HLA-DRA major histocompatibility complex, class II, DR alpha
3597981_st -1.7779721 1.957269 -13.205556 1.00e-06 0.0059623 5.058531 — —
47423221_st 1.2842290 2.008371 11.678282 1.50e-06 0.0059623 4.982022 — —
TC6_mann_hap4000074.hg.1 1.5788721 5.509731 11.618619 1.50e-06 0.0059623 4.954878 HLA-DRA major histocompatibility complex, class II, DR alpha
TC6_mcf_hap5000075.hg.1 1.5788721 5.509731 11.618619 1.50e-06 0.0059623 4.954878 HLA-DRA major histocompatibility complex, class II, DR alpha
TC6_qbl_hap6000081.hg.1 1.5788721 5.509731 11.618619 1.50e-06 0.0059623 4.954878 HLA-DRA major histocompatibility complex, class II, DR alpha
TC6_ssto_hap7000076.hg.1 1.5788721 5.509731 11.618619 1.50e-06 0.0059623 4.954878 HLA-DRA major histocompatibility complex, class II, DR alpha
47420668_st -1.2842484 2.570708 -10.990091 2.40e-06 0.0083630 4.654413 — —
TC06001902.hg.1 1.6047891 5.493022 10.944078 2.50e-06 0.0083630 4.631331 ELOVL4 ELOVL fatty acid elongase 4
TC09000561.hg.1 1.3587189 3.064936 10.625651 3.10e-06 0.0100387 4.467307 MIR4668 microRNA 4668
TC02002218.hg.1 2.1814536 5.332032 12.295575 2.50e-06 0.0083630 4.431663 IL1A interleukin 1, alpha
TC0X001832.hg.1 1.3065124 3.106171 10.491845 3.50e-06 0.0100549 4.396082 — —
TC08001992.hg.1 1.2430949 2.601038 10.465772 3.50e-06 0.0100549 4.382040 — —
TC0X001951.hg.1 1.4028290 5.097887 10.455654 3.60e-06 0.0100549 4.376576 — —
TC15002198.hg.1 1.3734181 3.049986 10.217507 4.30e-06 0.0116265 4.245616 TMOD2 tropomodulin 2 (neuronal)
47422683_st -1.2069796 2.322069 -10.022148 5.00e-06 0.0126628 4.134708 — —
TC11003477.hg.1 1.3242186 4.456241 10.015378 5.00e-06 0.0126628 4.130807 SAA4 serum amyloid A4, constitutive
TC12002256.hg.1 1.5074796 6.362556 9.951028 5.30e-06 0.0128704 4.093532 — —
2985881_st -1.2580780 4.564692 -9.832541 5.80e-06 0.0136855 4.023961 — —
TC03000407.hg.1 1.2284360 4.285724 9.635331 6.80e-06 0.0155547 3.905420 KBTBD8 kelch repeat and BTB (POZ) domain containing 8
TC01003017.hg.1 1.3864164 6.241617 9.514262 7.60e-06 0.0164101 3.830904 TSPAN2 tetraspanin 2
TC12002778.hg.1 1.2756243 5.401293 9.495111 7.70e-06 0.0164101 3.818993 MGP matrix Gla protein
TC05000159.hg.1 1.3522513 6.009540 9.384238 8.40e-06 0.0169245 3.749362 IL7R interleukin 7 receptor
TC02004676.hg.1 1.2157641 3.686998 9.369066 8.50e-06 0.0169245 3.739743 SCN9A sodium channel, voltage-gated, type IX, alpha subunit
TC11002081.hg.1 1.1960821 4.538888 9.301507 9.00e-06 0.0169245 3.696647 PGM2L1 phosphoglucomutase 2-like 1
TC03002268.hg.1 1.2169468 4.237264 9.284619 9.20e-06 0.0169245 3.685805 TGFBR2 transforming growth factor, beta receptor II (70/80kDa)
TC12000189.hg.1 1.3289915 6.172095 9.263642 9.30e-06 0.0169245 3.672301 EMP1 epithelial membrane protein 1
TC08001148.hg.1 1.1549631 3.809360 9.259949 9.40e-06 0.0169245 3.669919 RPS20P22 ribosomal protein S20 pseudogene 22
TC05000084.hg.1 1.2531189 3.872873 9.131389 1.04e-05 0.0184198 3.586171 FAM105A family with sequence similarity 105, member A
47421342_st -1.1627233 2.507763 -9.085310 1.09e-05 0.0186990 3.555759 — —
47424037_st -1.0015761 2.120592 -8.916886 1.26e-05 0.0211374 3.442782 — —
2839532_st -1.3517704 6.952548 -8.826167 1.36e-05 0.0219484 3.380728 — —
TC08001361.hg.1 1.1949998 5.260314 8.808719 1.38e-05 0.0219484 3.368696 PAG1 phosphoprotein membrane anchor with glycosphingolipid microdomains 1
TC0X001577.hg.1 1.0934401 3.869794 8.739513 1.47e-05 0.0225791 3.320656 — —
TC08001112.hg.1 1.1217793 4.654709 8.677915 1.56e-05 0.0230287 3.277471 TEX15 testis expressed 15
TC0X001575.hg.1 1.1463443 5.172118 8.669945 1.57e-05 0.0230287 3.271854 — —
TC20000642.hg.1 1.2142711 3.983362 8.627762 1.63e-05 0.0234314 3.242011 FLRT3 fibronectin leucine rich transmembrane protein 3
TC08002280.hg.1 1.8270632 3.707449 10.269895 1.40e-05 0.0219484 3.225239 — —
TC05002579.hg.1 1.1668672 3.979395 8.565839 1.72e-05 0.0242854 3.197855 — —
TC11003478.hg.1 1.2346532 6.467708 8.409094 1.99e-05 0.0270476 3.084211 SAA2 serum amyloid A2
TC6_mcf_hap5000068.hg.1 -1.0443836 3.898554 -8.354081 2.09e-05 0.0273159 3.043679 C4B complement component 4B (Childo blood group)
TC04001347.hg.1 1.2559545 5.049703 8.353480 2.09e-05 0.0273159 3.043234 HPSE heparanase
47423273_st -1.3644658 2.381801 -9.156331 1.99e-05 0.0270476 3.031178 — —
TC09002617.hg.1 1.2169156 4.710199 8.279808 2.24e-05 0.0286668 2.988415 RASEF RAS and EF-hand domain containing
47420572_st -1.0787343 2.587622 -8.261882 2.28e-05 0.0286668 2.974983 — —
TC0X000544.hg.1 1.0272285 3.411397 8.243492 2.32e-05 0.0286668 2.961165 — —
TC0X001610.hg.1 1.0933317 2.870326 8.181863 2.45e-05 0.0298495 2.914574 REPS2 RALBP1 associated Eps domain containing 2
TC08001988.hg.1 1.1400060 4.296443 8.130047 2.58e-05 0.0308136 2.875061 PKIA protein kinase (cAMP-dependent, catalytic) inhibitor alpha
TC10002196.hg.1 1.1814331 5.694986 8.086857 2.69e-05 0.0315661 2.841886 FAS Fas cell surface death receptor
TC06002615.hg.1 1.4628795 4.877372 8.709496 2.75e-05 0.0318377 2.788936 SOX4 SRY (sex determining region Y)-box 4
47423203_st 1.1279281 3.040097 8.045277 2.93e-05 0.0332918 2.769242 — —
3371087_st -1.0056388 3.956589 -7.980006 2.97e-05 0.0332918 2.758866 — —
TC09001842.hg.1 1.0175019 4.569505 7.929773 3.12e-05 0.0339668 2.719365 — —
TC12001226.hg.1 1.0666160 4.421865 7.867470 3.32e-05 0.0339668 2.669949 STYK1 serine/threonine/tyrosine kinase 1
TC06001401.hg.1 1.1127446 6.078436 7.851139 3.37e-05 0.0339668 2.656918 ZNF184 zinc finger protein 184
TC06002062.hg.1 1.1287863 4.909093 7.827256 3.45e-05 0.0339668 2.637802 MAN1A1 mannosidase, alpha, class 1A, member 1
47421855_st -0.8767796 2.157655 -7.820430 3.47e-05 0.0339668 2.632326 — —
TC6_mcf_hap5000213.hg.1 1.2231736 6.662812 7.798362 3.55e-05 0.0339668 2.614582 CFB complement factor B
TC6_ssto_hap7000199.hg.1 1.2231736 6.662812 7.798362 3.55e-05 0.0339668 2.614582 CFB complement factor B
TC6_dbb_hap3000222.hg.1 1.2209378 6.663288 7.789521 3.58e-05 0.0339668 2.607456 CFB complement factor B
TC12001755.hg.1 1.2113055 6.297977 7.787771 3.58e-05 0.0339668 2.606044 CSRP2 cysteine and glycine-rich protein 2
TC01001624.hg.1 1.0037725 4.219329 7.783327 3.60e-05 0.0339668 2.602458 RGS2 regulator of G-protein signaling 2
TC6_cox_hap2000245.hg.1 1.2208918 6.663264 7.781030 3.61e-05 0.0339668 2.600604 CFB complement factor B
TC6_qbl_hap6000225.hg.1 1.2176023 6.662276 7.769288 3.65e-05 0.0339668 2.591113 CFB complement factor B
TC03001888.hg.1 1.2857128 5.943545 7.854371 3.69e-05 0.0339668 2.578855 TM4SF1 transmembrane 4 L six family member 1
TC01002163.hg.1 1.3323629 7.444205 7.749341 3.72e-05 0.0339668 2.574950 TNFRSF9 tumor necrosis factor receptor superfamily, member 9
TC6_mann_hap4000196.hg.1 1.1694408 6.674467 7.736564 3.77e-05 0.0339668 2.564572 CFB complement factor B
TC22000762.hg.1 -0.9939818 4.466419 -7.726509 3.80e-05 0.0339668 2.556390 — —
TC17000720.hg.1 1.2075208 5.445204 7.732132 3.96e-05 0.0345529 2.522498 GDPD1 glycerophosphodiester phosphodiesterase domain containing 1
TC12000182.hg.1 1.1515985 6.484116 7.683819 3.97e-05 0.0345529 2.521512 GPRC5A G protein-coupled receptor, class C, group 5, member A
TC02003411.hg.1 1.0790584 4.877519 7.647460 4.11e-05 0.0353829 2.491626 — —
TC06004077.hg.1 1.1799352 6.560729 7.621929 4.22e-05 0.0358545 2.470540 CFB complement factor B
TC08001820.hg.1 1.1665489 6.004030 7.482157 4.85e-05 0.0403118 2.353631 NAT1 N-acetyltransferase 1 (arylamine N-acetyltransferase)
47423846_st -0.7673152 1.464740 -7.481196 4.86e-05 0.0403118 2.352818 — —
47422398_st -0.8970998 2.769027 -7.419277 5.17e-05 0.0418053 2.300212 — —
TC19002417.hg.1 0.9452021 4.141262 7.412380 5.21e-05 0.0418053 2.294321 — —
TC11002162.hg.1 0.9311288 4.022142 7.394992 5.30e-05 0.0418053 2.279443 SYTL2 synaptotagmin-like 2
TC02001140.hg.1 0.9272994 2.735180 7.389322 5.33e-05 0.0418053 2.274583 LOC101927482 uncharacterized LOC101927482
TC03003320.hg.1 1.1437088 6.573990 7.388926 5.34e-05 0.0418053 2.274243 NFKBIZ nuclear factor of kappa light polypeptide gene enhancer in B-cells inhibitor, zeta
TC02002747.hg.1 1.0629501 5.536560 7.327208 5.68e-05 0.0440371 2.221061 FN1 fibronectin 1
TC08002279.hg.1 1.2222416 5.172348 7.558510 5.84e-05 0.0440836 2.197066 — —
TC01005829.hg.1 0.9782791 5.381701 7.291486 5.89e-05 0.0440836 2.190048 FLG filaggrin
TC12001282.hg.1 1.0452450 6.190054 7.287076 5.92e-05 0.0440836 2.186208 EPS8 epidermal growth factor receptor pathway substrate 8
TC09002903.hg.1 1.1708804 6.395968 7.284305 5.94e-05 0.0440836 2.183794 TMEFF1 transmembrane protein with EGF-like and two follistatin-like domains 1
TC18000683.hg.1 0.9563513 3.378828 7.270776 6.02e-05 0.0442367 2.171991 — —
TC16000258.hg.1 -1.1259825 8.965091 -7.209801 6.41e-05 0.0463319 2.118491 PLK1 polo-like kinase 1
TC21000944.hg.1 0.8969179 3.698593 7.206107 6.44e-05 0.0463319 2.115233 — —
TC20001333.hg.1 -1.0191143 6.070298 -7.187519 6.56e-05 0.0465727 2.098814 — —
TC12003271.hg.1 0.9474799 3.527480 7.181677 6.60e-05 0.0465727 2.093644 KLRC3 killer cell lectin-like receptor subfamily C, member 3
TC14001730.hg.1 0.9188216 4.656683 7.162118 6.74e-05 0.0470606 2.076303 AKAP5 A kinase (PRKA) anchor protein 5
TC04001286.hg.1 1.0482472 5.317507 7.132672 6.95e-05 0.0480543 2.050095 CXCL2 chemokine (C-X-C motif) ligand 2
PSR17013959.hg.1 -0.9435930 5.145416 -7.093690 7.24e-05 0.0495719 2.015217 — —
TC04001698.hg.1 1.0195227 4.247191 7.066124 7.45e-05 0.0499656 1.990428 NPY1R neuropeptide Y receptor Y1
47420654_st -1.2099436 2.119399 -7.935327 7.43e-05 0.0499656 1.988252 — —
47419844_st -0.7441258 1.453461 -7.058879 7.51e-05 0.0499656 1.983896 — —

This table can be regarded as a final output of the pairwise differential expression pipeline.

9.split this into up- and down-regulated gene tables and lists for analysis and pathway mapping

up <- diff[diff$logFC > 0,]
write.csv(up,"pairwise_test_results/tGFR.cGFR/up.tGFR.cGFR",row.names = F,quote = F)
write.csv(up$name,"pairwise_test_results/tGFR.cGFR/up.tGFR.cGFR.gene.list",row.names = F,quote = F)
down <- diff[diff$logFC < 0,]
write.csv(down,"pairwise_test_results/tGFR.cGFR/down.tGFR.cGFR",row.names = F,quote = F)
write.csv(down$name,"pairwise_test_results/tGFR.cGFR/down.tGFR.cGFR.gene.list",row.names = F,quote = F)

10. Plotting intensity and variance, highlighting significant differences (plot shown at start of document also)

png("pairwise_test_results/tGFR.cGFR/tGFR.cGFR.MD.png")
plotMD(fit,status = results[,3], legend = F, hl.col = c("red","blue"), main = "tGFR vs. cGFR", xlim=c(0,14))
abline(0,0,col="red",lty=3)
dev.off()
## png 
##   2
plotMD(fit,status = results[,3], legend = F, hl.col = c("red","blue"), main = "tGFR vs. cGFR", xlim=c(0,14))
abline(0,0,col="red",lty=3)

We now have full tables for differentially expressed genes. We can also generate tables of up- and down-regulated pathways via the Reactome.db. These will be handled in a separate script.

11. Linking to GO enrichment analysis: step 1, gather a set of background genes to compare

This pipeline is adapted from https://bioconductor.org/packages/devel/workflows/vignettes/maEndToEnd/inst/doc/MA-Workflow.html (this required a good deal of modification to current pipeline)

First, we get a “control set” of genes, those which have similar expression levels but are not differentially expressed. this is handled by a package “genefilter”.

library(genefilter)
library(Biobase)
#make a list of diff genes
diff.probe <- subset(x = diff)$probe
background.genes <- genefilter::genefinder(eset, as.character(diff.probe), method = "manhattan", scale = "none")
background.genes <- sapply(background.genes, function(x) x$indices)
background.genes <- as.vector(background.genes)
back.genes <- featureNames(project.bgcorrect.norm.avg)[background.genes]

We now remove any differentially expressed genes from the background gene set and verify that it is correct: the diff set and the background set should have nothing in common.

back.genes <- setdiff(back.genes, diff.probe)
intersect(back.genes, diff.probe)
## character(0)

12. Examine mean intensity of the diff set, background set and all probes

library(cancerTiming)
multidensity(list(
    all = fit.tGFR.cGFR[,"AveExpr"],
    fore = diff[,"AveExpr"],  
    back = fit.tGFR.cGFR[rownames(fit.tGFR.cGFR) %in% back.genes, "AveExpr"]),
  col = c("#e46981", "#ae7ee2", "#a7ad4a"),
  xlab = "mean expression",
  main = "cGFR vs cGF: Differentially Expressed, Background, and All gene expression levels"
)

13. GO enrichment analysis

Is the differential expression set enriched in particular functional categories? This is why we have created a background gene set, so that we know the “typical” functional distribution of genes expressed at similar levels. We will connect the probe names to GO categories and run a few tests for pathway enrichment.

library(topGO)

gene.IDs <- rownames(eset)
in.universe <- gene.IDs %in% c(diff.probe,back.genes)
in.selection <- gene.IDs %in% diff.probe

all.genes <- in.selection[in.universe]
all.genes <- factor(as.integer(in.selection[in.universe]))
names(all.genes) <- gene.IDs[in.universe]

#this builds a database with all GO information related to our two test sets AND whether the gene is differentially expressed
top_GO_data <- new("topGOdata", ontology = "BP", allGenes = all.genes,
 nodeSize = 10, annot = annFUN.db, affyLib = "hta20transcriptcluster.db")

result_top_GO_elim <- 
  runTest(top_GO_data, algorithm = "elim", statistic = "Fisher")
result_top_GO_classic <- 
  runTest(top_GO_data, algorithm = "classic", statistic = "Fisher")
library(stringr)
res_top_GO <- GenTable(top_GO_data, Fisher.elim = result_top_GO_elim,
        Fisher.classic = result_top_GO_classic,
        orderBy = "Fisher.elim" , topNodes = 100)

genes_top_GO <- printGenes(top_GO_data, whichTerms = res_top_GO$GO.ID,
    chip = "hta20transcriptcluster.db", geneCutOff = 1000)

res_top_GO$sig_genes <- sapply(genes_top_GO, function(x){
                str_c(paste0(x[x$'raw p-value' == 2, "Symbol.id"],";"), 
                      collapse = "")
    })
kable(res_top_GO) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")
GO.ID Term Annotated Significant Expected Rank in Fisher.classic Fisher.elim Fisher.classic sig_genes
GO:0050852 T cell receptor signaling pathway 15 11 2.48 3 8.9e-07 8.9e-07 NFKBIZ;HLA-DRA;PAG1;PAK3;PAK3;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;
GO:0050890 cognition 10 8 1.65 9 1.3e-05 1.3e-05 HLA-DRA;TMOD2;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;
GO:0002673 regulation of acute inflammatory respons… 10 7 1.65 16 0.00020 0.00020 CFB;CFB;CFB;CFB;CFB;CFB;CFB;
GO:0019882 antigen processing and presentation 10 7 1.65 17 0.00020 0.00020 HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;
GO:0071346 cellular response to interferon-gamma 11 7 1.82 23 0.00047 0.00047 HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;
GO:0070613 regulation of protein processing 11 7 1.82 24 0.00047 0.00047 CFB;CFB;CFB;CFB;CFB;CFB;CFB;
GO:0034622 cellular protein-containing complex asse… 26 11 4.30 27 0.00109 0.00109 HLA-DRA;PAK3;PAK3;EPS8;TMOD2;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;
GO:0019221 cytokine-mediated signaling pathway 32 12 5.29 28 0.00226 0.00226 TNFRSF9;IL1A;FN1;CXCL2;IL7R;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;
GO:0002697 regulation of immune effector process 21 9 3.47 33 0.00300 0.00300 NFKBIZ;IL7R;CFB;CFB;CFB;CFB;CFB;CFB;CFB;
GO:0002253 activation of immune response 31 18 5.12 1 0.00312 3.2e-08 NFKBIZ;HLA-DRA;CFB;PAG1;PAK3;PAK3;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;
GO:0002526 acute inflammatory response 16 11 2.64 6 0.00549 2.5e-06 IL1A;FN1;CFB;SAA4;SAA2;CFB;CFB;CFB;CFB;CFB;CFB;
GO:0006959 humoral immune response 19 8 3.14 38 0.00603 0.00603 CXCL2;CFB;CFB;CFB;CFB;CFB;CFB;CFB;
GO:0044403 symbiont process 20 8 3.31 41 0.00874 0.00874 FN1;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;
GO:0090066 regulation of anatomical structure size 18 7 2.98 47 0.01711 0.01711 RGS2;FN1;IL7R;PAK3;PAK3;EPS8;TMOD2;
GO:0045087 innate immune response 40 16 6.61 13 0.01832 0.00013 HLA-DRA;CFB;PAK3;PAK3;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;
GO:0016032 viral process 19 7 3.14 50 0.02358 0.02358 HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;
GO:0032535 regulation of cellular component size 15 6 2.48 52 0.02374 0.02374 FN1;IL7R;PAK3;PAK3;EPS8;TMOD2;
GO:1903037 regulation of leukocyte cell-cell adhesi… 18 6 2.98 57 0.05836 0.05836 TGFBR2;NFKBIZ;IL7R;PAG1;PAK3;PAK3;
GO:0050863 regulation of T cell activation 18 6 2.98 58 0.05836 0.05836 TGFBR2;NFKBIZ;IL7R;PAG1;PAK3;PAK3;
GO:0007159 leukocyte cell-cell adhesion 18 6 2.98 59 0.05836 0.05836 TGFBR2;NFKBIZ;IL7R;PAG1;PAK3;PAK3;
GO:0033673 negative regulation of kinase activity 10 4 1.65 61 0.06531 0.06531 RGS2;PKIA;NA;PLK1;
GO:0110053 regulation of actin filament organizatio… 10 4 1.65 62 0.06531 0.06531 PAK3;PAK3;EPS8;TMOD2;
GO:0051052 regulation of DNA metabolic process 10 4 1.65 63 0.06531 0.06531 IL7R;TEX15;PAK3;PAK3;
GO:0051348 negative regulation of transferase activ… 10 4 1.65 64 0.06531 0.06531 RGS2;PKIA;NA;PLK1;
GO:0006469 negative regulation of protein kinase ac… 10 4 1.65 65 0.06531 0.06531 RGS2;PKIA;NA;PLK1;
GO:0042110 T cell activation 23 7 3.80 66 0.06587 0.06587 TGFBR2;NFKBIZ;IL7R;SOX4;PAG1;PAK3;PAK3;
GO:0051249 regulation of lymphocyte activation 19 6 3.14 67 0.07441 0.07441 TGFBR2;NFKBIZ;IL7R;PAG1;PAK3;PAK3;
GO:0030155 regulation of cell adhesion 29 8 4.79 69 0.08458 0.08458 FN1;TGFBR2;NFKBIZ;IL7R;PAG1;ACER2;PAK3;PAK3;
GO:0022407 regulation of cell-cell adhesion 20 6 3.31 71 0.09277 0.09277 TGFBR2;NFKBIZ;IL7R;PAG1;PAK3;PAK3;
GO:0044087 regulation of cellular component biogene… 20 6 3.31 72 0.09277 0.09277 PAK3;PAK3;EPS8;TMOD2;PLK1;FLRT3;
GO:0045785 positive regulation of cell adhesion 20 6 3.31 73 0.09277 0.09277 FN1;TGFBR2;NFKBIZ;IL7R;PAK3;PAK3;
GO:1903039 positive regulation of leukocyte cell-ce… 16 5 2.64 76 0.10561 0.10561 TGFBR2;NFKBIZ;IL7R;PAK3;PAK3;
GO:0051251 positive regulation of lymphocyte activa… 16 5 2.64 77 0.10561 0.10561 TGFBR2;NFKBIZ;IL7R;PAK3;PAK3;
GO:0050870 positive regulation of T cell activation 16 5 2.64 78 0.10561 0.10561 TGFBR2;NFKBIZ;IL7R;PAK3;PAK3;
GO:0098609 cell-cell adhesion 26 7 4.30 79 0.11687 0.11687 TGFBR2;NFKBIZ;IL7R;PAG1;PAK3;PAK3;FLRT3;
GO:0010976 positive regulation of neuron projection… 12 4 1.98 80 0.11885 0.11885 RGS2;FN1;PAK3;PAK3;
GO:0007015 actin filament organization 12 4 1.98 81 0.11885 0.11885 PAK3;PAK3;EPS8;TMOD2;
GO:0032956 regulation of actin cytoskeleton organiz… 12 4 1.98 82 0.11885 0.11885 PAK3;PAK3;EPS8;TMOD2;
GO:1902903 regulation of supramolecular fiber organ… 12 4 1.98 83 0.11885 0.11885 PAK3;PAK3;EPS8;TMOD2;
GO:0032970 regulation of actin filament-based proce… 12 4 1.98 84 0.11885 0.11885 PAK3;PAK3;EPS8;TMOD2;
GO:0043604 amide biosynthetic process 12 4 1.98 85 0.11885 0.11885 RGS2;KBTBD8;ELOVL4;SOX4;
GO:0048588 developmental cell growth 12 4 1.98 86 0.11885 0.11885 RGS2;FN1;TGFBR2;FLRT3;
GO:0022610 biological adhesion 42 10 6.94 87 0.13098 0.13098 FN1;TGFBR2;NFKBIZ;HPSE;IL7R;PAG1;ACER2;PAK3;PAK3;FLRT3;
GO:0007155 cell adhesion 42 10 6.94 88 0.13098 0.13098 FN1;TGFBR2;NFKBIZ;HPSE;IL7R;PAG1;ACER2;PAK3;PAK3;FLRT3;
GO:0031346 positive regulation of cell projection o… 17 5 2.81 89 0.13103 0.13103 RGS2;FN1;PAK3;PAK3;EPS8;
GO:0022409 positive regulation of cell-cell adhesio… 17 5 2.81 90 0.13103 0.13103 TGFBR2;NFKBIZ;IL7R;PAK3;PAK3;
GO:0043603 cellular amide metabolic process 17 5 2.81 91 0.13103 0.13103 RGS2;KBTBD8;ELOVL4;SOX4;ACER2;
GO:0051493 regulation of cytoskeleton organization 17 5 2.81 92 0.13103 0.13103 PAK3;PAK3;EPS8;TMOD2;PLK1;
GO:0002694 regulation of leukocyte activation 22 6 3.64 93 0.13619 0.13619 TGFBR2;NFKBIZ;IL7R;PAG1;PAK3;PAK3;
GO:0050865 regulation of cell activation 22 6 3.64 94 0.13619 0.13619 TGFBR2;NFKBIZ;IL7R;PAG1;PAK3;PAK3;
GO:0007610 behavior 13 4 2.15 96 0.15102 0.15102 SCN9A;NPY1R;EPS8;TMOD2;
GO:0006952 defense response 63 25 10.41 2 0.16251 7.3e-07 TSPAN2;IL1A;FN1;SCN9A;NFKBIZ;CXCL2;HLA-DRA;CFB;PAK3;PAK3;SAA4;SAA2;KLRC3;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;
GO:0070887 cellular response to chemical stimulus 89 23 14.71 39 0.17924 0.00672 TNFRSF9;IL1A;FN1;TGFBR2;CXCL2;IL7R;HLA-DRA;SOX4;NAT1;ACER2;PAK3;PAK3;FAS;SAA4;SAA2;EPS8;FLRT3;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;HLA-DRA;
GO:0046649 lymphocyte activation 29 7 4.79 98 0.18388 0.18388 TGFBR2;NFKBIZ;IL7R;SOX4;PAG1;PAK3;PAK3;
GO:0045862 positive regulation of proteolysis 14 4 2.31 100 0.18614 0.18614 FN1;ACER2;FAS;PLK1;
GO:0048638 regulation of developmental growth 14 4 2.31 101 0.18614 0.18614 RGS2;FN1;TGFBR2;NPY1R;
GO:0045666 positive regulation of neuron differenti… 14 4 2.31 102 0.18614 0.18614 RGS2;FN1;PAK3;PAK3;
GO:0051247 positive regulation of protein metabolic… 45 10 7.44 103 0.18630 0.18630 IL1A;FN1;TGFBR2;SOX4;ACER2;PAK3;PAK3;FAS;NA;PLK1;
GO:0002696 positive regulation of leukocyte activat… 19 5 3.14 104 0.18918 0.18918 TGFBR2;NFKBIZ;IL7R;PAK3;PAK3;
GO:0050867 positive regulation of cell activation 19 5 3.14 105 0.18918 0.18918 TGFBR2;NFKBIZ;IL7R;PAK3;PAK3;
GO:0002376 immune system process 98 28 16.20 20 0.19241 0.00028 TSPAN2;IL1A;FN1;TGFBR2;NFKBIZ;CXCL2;HPSE;IL7R;HLA-DRA;SOX4;CFB;PAG1;PAK3;PAK3;FAS;EPS8;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;HLA-DRA;CFB;
GO:1901566 organonitrogen compound biosynthetic pro… 35 8 5.79 107 0.20190 0.20190 RGS2;IL1A;KBTBD8;ELOVL4;MAN1A1;SOX4;ACER2;PGM2L1;
GO:0006793 phosphorus metabolic process 74 15 12.23 108 0.21124 0.21124 RGS2;FN1;TGFBR2;ELOVL4;PKIA;PAK3;PAK3;FAS;PGM2L1;SYTL2;SYTL2;NA;STYK1;PLK1;GDPD1;
GO:0006796 phosphate-containing compound metabolic … 74 15 12.23 109 0.21124 0.21124 RGS2;FN1;TGFBR2;ELOVL4;PKIA;PAK3;PAK3;FAS;PGM2L1;SYTL2;SYTL2;NA;STYK1;PLK1;GDPD1;
GO:0007169 transmembrane receptor protein tyrosine … 25 6 4.13 110 0.21601 0.21601 PAK3;REPS2;PAK3;NA;EPS8;FLRT3;
GO:0007188 adenylate cyclase-modulating G protein-c… 10 3 1.65 111 0.21902 0.21902 RGS2;NPY1R;AKAP5;
GO:0097305 response to alcohol 10 3 1.65 112 0.21902 0.21902 RGS2;TGFBR2;EPS8;
GO:0043043 peptide biosynthetic process 10 3 1.65 113 0.21902 0.21902 RGS2;KBTBD8;SOX4;
GO:0050808 synapse organization 10 3 1.65 114 0.21902 0.21902 PAK3;PAK3;FLRT3;
GO:0031589 cell-substrate adhesion 10 3 1.65 115 0.21902 0.21902 FN1;HPSE;ACER2;
GO:0006412 translation 10 3 1.65 116 0.21902 0.21902 RGS2;KBTBD8;SOX4;
GO:0006417 regulation of translation 10 3 1.65 117 0.21902 0.21902 RGS2;KBTBD8;SOX4;
GO:0010952 positive regulation of peptidase activit… 10 3 1.65 118 0.21902 0.21902 FN1;ACER2;FAS;
GO:0007187 G protein-coupled receptor signaling pat… 10 3 1.65 119 0.21902 0.21902 RGS2;NPY1R;AKAP5;
GO:0071900 regulation of protein serine/threonine k… 20 5 3.31 120 0.22132 0.22132 RGS2;PKIA;PAK3;PAK3;PLK1;
GO:0048589 developmental growth 20 5 3.31 121 0.22132 0.22132 RGS2;FN1;TGFBR2;NPY1R;FLRT3;
GO:0000280 nuclear division 15 4 2.48 122 0.22367 0.22367 IL1A;TEX15;EPS8;PLK1;
GO:0032147 activation of protein kinase activity 15 4 2.48 123 0.22367 0.22367 TGFBR2;PAK3;PAK3;NA;
GO:0048598 embryonic morphogenesis 15 4 2.48 124 0.22367 0.22367 FN1;TGFBR2;SOX4;FLRT3;
GO:0045859 regulation of protein kinase activity 31 7 5.12 128 0.23602 0.23602 RGS2;TGFBR2;PKIA;PAK3;PAK3;NA;PLK1;
GO:0043549 regulation of kinase activity 31 7 5.12 129 0.23602 0.23602 RGS2;TGFBR2;PKIA;PAK3;PAK3;NA;PLK1;
GO:0051174 regulation of phosphorus metabolic proce… 53 11 8.76 130 0.23854 0.23854 RGS2;FN1;TGFBR2;PKIA;PAK3;PAK3;FAS;SYTL2;SYTL2;NA;PLK1;
GO:0019220 regulation of phosphate metabolic proces… 53 11 8.76 131 0.23854 0.23854 RGS2;FN1;TGFBR2;PKIA;PAK3;PAK3;FAS;SYTL2;SYTL2;NA;PLK1;
GO:0043086 negative regulation of catalytic activit… 26 6 4.30 133 0.24575 0.24575 RGS2;PKIA;SYTL2;SYTL2;NA;PLK1;
GO:0045936 negative regulation of phosphate metabol… 26 6 4.30 134 0.24575 0.24575 RGS2;PKIA;SYTL2;SYTL2;NA;PLK1;
GO:0010563 negative regulation of phosphorus metabo… 26 6 4.30 135 0.24575 0.24575 RGS2;PKIA;SYTL2;SYTL2;NA;PLK1;
GO:0051246 regulation of protein metabolic process 82 20 13.56 53 0.24907 0.02535 RGS2;IL1A;FN1;KBTBD8;TGFBR2;SOX4;CFB;PKIA;ACER2;PAK3;PAK3;FAS;NA;PLK1;CFB;CFB;CFB;CFB;CFB;CFB;
GO:0051962 positive regulation of nervous system de… 21 5 3.47 136 0.25504 0.25504 RGS2;FN1;PAK3;PAK3;FLRT3;
GO:0030217 T cell differentiation 16 4 2.64 137 0.26302 0.26302 TGFBR2;NFKBIZ;IL7R;SOX4;
GO:0048285 organelle fission 16 4 2.64 138 0.26302 0.26302 IL1A;TEX15;EPS8;PLK1;
GO:0051130 positive regulation of cellular componen… 32 7 5.29 139 0.26385 0.26385 RGS2;IL1A;FN1;PAK3;PAK3;EPS8;FLRT3;
GO:0032270 positive regulation of cellular protein … 43 9 7.11 140 0.26400 0.26400 FN1;TGFBR2;SOX4;ACER2;PAK3;PAK3;FAS;NA;PLK1;
GO:0045621 positive regulation of lymphocyte differ… 11 3 1.82 141 0.26735 0.26735 TGFBR2;NFKBIZ;IL7R;
GO:0045619 regulation of lymphocyte differentiation 11 3 1.82 142 0.26735 0.26735 TGFBR2;NFKBIZ;IL7R;
GO:0034248 regulation of cellular amide metabolic p… 11 3 1.82 143 0.26735 0.26735 RGS2;KBTBD8;SOX4;
GO:0045580 regulation of T cell differentiation 11 3 1.82 144 0.26735 0.26735 TGFBR2;NFKBIZ;IL7R;
GO:0045582 positive regulation of T cell differenti… 11 3 1.82 145 0.26735 0.26735 TGFBR2;NFKBIZ;IL7R;
GO:0042327 positive regulation of phosphorylation 27 6 4.46 146 0.27667 0.27667 TGFBR2;PAK3;PAK3;FAS;NA;PLK1;
GO:0001934 positive regulation of protein phosphory… 27 6 4.46 147 0.27667 0.27667 TGFBR2;PAK3;PAK3;FAS;NA;PLK1;
GO:0033043 regulation of organelle organization 27 6 4.46 148 0.27667 0.27667 IL1A;PAK3;PAK3;EPS8;TMOD2;PLK1; 
write.csv(res_top_GO, "pairwise_test_results/tGFR.cGFR/tGFR.cGFR.GO.enrichment.csv")

This is a ranking of which GO terms are enriched according to two different testing strategies, elim and Fisher.classic. these are, roughly, probability scores and can be broadly interpreted the same way as p-scores.

14. Visualizing the results:

nodes <- showSigOfNodes(top_GO_data, score(result_top_GO_elim), firstSigNodes = 3,
               useInfo = 'def')

png("pairwise_test_results/tGFR.cGFR/tGFR.tGFC.GO.nodes.png")
nodes
## $dag
## A graphNEL graph with directed edges
## Number of Nodes = 41 
## Number of Edges = 71 
## 
## $complete.dag
## [1] "A graph with 41 nodes."
dev.off()
## png 
##   2
nodes
## $dag
## A graphNEL graph with directed edges
## Number of Nodes = 41 
## Number of Edges = 71 
## 
## $complete.dag
## [1] "A graph with 41 nodes."

15. Pathway enrichment analysis via Reactome

Reactome maps

library(reactome.db)
library(DOSE)
library(ReactomePA)
entrez_ids <- mapIds(hta20transcriptcluster.db, 
      keys = rownames(eset), 
      keytype = "PROBEID",
      column = "ENTREZID",
      multiVals = "first")

reactome_enrich <- enrichPathway(gene = entrez_ids[diff.probe], 
                                universe = entrez_ids[c(diff.probe, 
                                                        back.genes)],
                                organism = "human",
                                pvalueCutoff = 1,
                                qvalueCutoff = 1, 
                                readable = TRUE)

reactome_enrich@result$Description <- paste0(str_sub(
                                    reactome_enrich@result$Description, 1, 20),
                                    "...")

write.csv(as.data.frame(reactome_enrich), "pairwise_test_results/tGFR.cGFR/tGFR.cGFR.reactome.enrichment.csv")

kable(reactome_enrich) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")
ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
R-HSA-9006934 R-HSA-9006934 Signaling by Recepto… 4/30 10/197 0.0478002 0.6726887 0.6726887 PAK3/PAG1/FLRT3/FN1 4
R-HSA-422475 R-HSA-422475 Axon guidance… 4/30 11/197 0.0670860 0.6726887 0.6726887 PAK3/SCN9A/FLRT3/AKAP5 4
R-HSA-1280215 R-HSA-1280215 Cytokine Signaling i… 6/30 21/197 0.0760266 0.6726887 0.6726887 HLA-DRA/IL1A/IL7R/TNFRSF9/FN1/CXCL2 6
R-HSA-1266738 R-HSA-1266738 Developmental Biolog… 5/30 17/197 0.0942811 0.6726887 0.6726887 PAK3/SCN9A/FLRT3/FLG/AKAP5 5
R-HSA-168256 R-HSA-168256 Immune System… 11/30 52/197 0.1240284 0.6726887 0.6726887 HLA-DRA/IL1A/PAK3/KBTBD8/IL7R/PAG1/HPSE/CFB/TNFRSF9/FN1/CXCL2 11
R-HSA-449147 R-HSA-449147 Signaling by Interle… 4/30 14/197 0.1452181 0.6726887 0.6726887 IL1A/IL7R/FN1/CXCL2 4
R-HSA-1430728 R-HSA-1430728 Metabolism… 7/30 31/197 0.1650451 0.6726887 0.6726887 ELOVL4/PGM2L1/HPSE/ACER2/GDPD1/NAT1/AKAP5 7
R-HSA-199991 R-HSA-199991 Membrane Trafficking… 3/30 10/197 0.1813722 0.6726887 0.6726887 IL7R/REPS2/MAN1A1 3
R-HSA-556833 R-HSA-556833 Metabolism of lipids… 3/30 11/197 0.2242296 0.6726887 0.6726887 ELOVL4/ACER2/GDPD1 3
R-HSA-5653656 R-HSA-5653656 Vesicle-mediated tra… 3/30 11/197 0.2242296 0.6726887 0.6726887 IL7R/REPS2/MAN1A1 3
R-HSA-1280218 R-HSA-1280218 Adaptive Immune Syst… 4/30 19/197 0.3229203 0.8806919 0.8806919 HLA-DRA/PAK3/KBTBD8/PAG1 4
R-HSA-162582 R-HSA-162582 Signal Transduction… 11/30 65/197 0.3939942 0.9712066 0.9712066 PAK3/TGFBR2/PAG1/FLRT3/FAS/SOX4/RGS2/FN1/PLK1/CXCL2/NPY1R 11
R-HSA-373076 R-HSA-373076 Class A/1 (Rhodopsin… 2/30 11/197 0.5216799 0.9712066 0.9712066 CXCL2/NPY1R 2
R-HSA-5683057 R-HSA-5683057 MAPK family signalin… 2/30 11/197 0.5216799 0.9712066 0.9712066 PAK3/FN1 2
R-HSA-194315 R-HSA-194315 Signaling by Rho GTP… 2/30 12/197 0.5722145 0.9712066 0.9712066 PAK3/PLK1 2
R-HSA-418594 R-HSA-418594 G alpha (i) signalli… 2/30 14/197 0.6619702 0.9712066 0.9712066 CXCL2/NPY1R 2
R-HSA-500792 R-HSA-500792 GPCR ligand binding… 2/30 14/197 0.6619702 0.9712066 0.9712066 CXCL2/NPY1R 2
R-HSA-597592 R-HSA-597592 Post-translational p… 4/30 29/197 0.6826186 0.9712066 0.9712066 KBTBD8/TGFBR2/MAN1A1/FN1 4
R-HSA-388396 R-HSA-388396 GPCR downstream sign… 3/30 23/197 0.7189955 0.9712066 0.9712066 RGS2/CXCL2/NPY1R 3
R-HSA-372790 R-HSA-372790 Signaling by GPCR… 3/30 24/197 0.7480967 0.9712066 0.9712066 RGS2/CXCL2/NPY1R 3
R-HSA-168249 R-HSA-168249 Innate Immune System… 3/30 25/197 0.7748887 0.9712066 0.9712066 PAK3/HPSE/CFB 3
R-HSA-1474244 R-HSA-1474244 Extracellular matrix… 1/30 10/197 0.8163330 0.9712066 0.9712066 FN1 1
R-HSA-1643685 R-HSA-1643685 Disease… 2/30 19/197 0.8231645 0.9712066 0.9712066 TGFBR2/FN1 2
R-HSA-392499 R-HSA-392499 Metabolism of protei… 4/30 35/197 0.8278317 0.9712066 0.9712066 KBTBD8/TGFBR2/MAN1A1/FN1 4
R-HSA-69278 R-HSA-69278 Cell Cycle, Mitotic… 1/30 11/197 0.8457983 0.9712066 0.9712066 PLK1 1
R-HSA-1640170 R-HSA-1640170 Cell Cycle… 1/30 12/197 0.8706695 0.9712066 0.9712066 PLK1 1
R-HSA-109582 R-HSA-109582 Hemostasis… 1/30 15/197 0.9241764 0.9712066 0.9712066 FN1 1
R-HSA-212436 R-HSA-212436 Generic Transcriptio… 2/30 28/197 0.9533076 0.9712066 0.9712066 FAS/ZNF184 2
R-HSA-73857 R-HSA-73857 RNA Polymerase II Tr… 2/30 31/197 0.9712066 0.9712066 0.9712066 FAS/ZNF184 2
R-HSA-74160 R-HSA-74160 Gene expression (Tra… 2/30 31/197 0.9712066 0.9712066 0.9712066 FAS/ZNF184 2

16. Visualizing reactome results

Barplots, showing reaction categories represented and their respective p-values:

png("pairwise_test_results/tGFR.cGFR/tGFR.cGFR.reactome.barplot.png")
barplot(reactome_enrich)
dev.off()
## png 
##   2
barplot(reactome_enrich)

png("pairwise_test_results/tGFR.cGFR/tGFR.cGFR.reactome.emapplot.png")
emapplot(reactome_enrich, showCategory = 10)
dev.off()
## png 
##   2
emapplot(reactome_enrich, showCategory = 10)

In this plot, color is p-value, size is proportional to the number of DE genes that mapped to the pathway, and the edge thickness is proportional to the number of shared genes between the units. This is useful for examining the interdependence of the pathways; which results should be taken more serioiusly, are independent? This particular plot shows the top ten most significant pathways, regardless of their p-values