Source: TBA

Introduction

Will conduct mitch analysis on reactome pathways for scRNA-seq data. Comparisons include naive vs bystander, bystander vs latent and latent vs productive. The focus is on apoptosis pathways and their member genes.

suppressPackageStartupMessages({
  library("mitch")
  library("gplots")
  library("kableExtra")
})

CORES=16

Loading processed data

The following comparisons were analyzed:

  1. latent_vs_bystander_markers.csv

  2. latent_vs_prod_markers.csv

  3. productive_vs_all_markers.csv

  4. productive_vs_mock_markers.csv

  5. latent_vs_mock_markers.csv

  6. Multi-dimensional (incl 1-5)

Note that in these comparisons, the 1st term is the “case” and the 2nd term is the “control”.

lat_v_by <- read.csv("latent_vs_bystander_markers.csv",row.names=1)

lat_v_prod <- read.csv("latent_vs_productive_markers.csv",row.names=1)

prod_v_all <- read.csv("productive_vs_all_markers.csv",row.names=1)

prod_v_mock <- read.csv("productive_vs_mock_markers.csv",row.names=1)

lat_v_mock <- read.csv("latent_vs_mock_markers.csv",row.names=1)

Gene sets

Anna, please check the spelling of the gene names.

gset <- gmt_import("ReactomePathways_2023-09-20.gmt")

apop <- c("BID","BCL2","BAD","BAK1","BAX","BBC3","BCL2L1","BCL2L2","BCL2L11","BCL2L13",
  "BCL2L14", "HRK", "MCL1", "BOK", "CYCS", "CCL5", "CCR5", "CDKN2A", "GHSR", "Gm14461",
  "MEF2C", "NOD2", "PLEKHO2", "PTEN", "SELENOS", "SIRT1", "ST6GAL1", "CASP3", "CASP7",
  "CASP8", "CASP9")

apop
##  [1] "BID"     "BCL2"    "BAD"     "BAK1"    "BAX"     "BBC3"    "BCL2L1" 
##  [8] "BCL2L2"  "BCL2L11" "BCL2L13" "BCL2L14" "HRK"     "MCL1"    "BOK"    
## [15] "CYCS"    "CCL5"    "CCR5"    "CDKN2A"  "GHSR"    "Gm14461" "MEF2C"  
## [22] "NOD2"    "PLEKHO2" "PTEN"    "SELENOS" "SIRT1"   "ST6GAL1" "CASP3"  
## [29] "CASP7"   "CASP8"   "CASP9"

Mitch analysis

1. Latent (case) vs bystander (ctrl).

lat_v_by

m1 <- mitch_import(lat_v_by,DEtype="seurat")
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 10939
## Note: no. genes in output = 10939
## Note: estimated proportion of input genes in output = 1
mres1 <- mitch_calc(x=m1,genesets=gset,cores=CORES,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
head(subset(mres1$enrichment_result,s.dist>0),20) %>%
  kbl(caption="pathways upregulated in latent compared to bystander") %>%
  kable_paper("hover",full_width=FALSE)
pathways upregulated in latent compared to bystander
set setSize pANOVA s.dist p.adjustANOVA
311 Erythropoietin activates Phosphoinositide-3-kinase (PI3K) 11 0.0006482 0.5939205 0.0032408
174 Cohesin Loading onto Chromatin 10 0.0015531 0.5779669 0.0072628
627 Mitotic Telophase/Cytokinesis 10 0.0015531 0.5779669 0.0072628
635 Myogenesis 15 0.0002546 0.5455145 0.0013508
911 Regulation of PTEN mRNA translation 11 0.0017694 0.5444563 0.0080710
655 NRAGE signals death through JNK 36 0.0000002 0.5046521 0.0000013
435 HDMs demethylate histones 21 0.0000683 0.5020717 0.0003945
178 Competing endogenous RNAs (ceRNAs) regulate PTEN translation 10 0.0060072 0.5018025 0.0221228
335 FOXO-mediated transcription of cell cycle genes 11 0.0039619 0.5017636 0.0157026
896 Regulation of CDH11 Expression and Function 15 0.0008500 0.4975711 0.0040933
1124 Synthesis of PIPs at the late endosome membrane 10 0.0069741 0.4927990 0.0248393
908 Regulation of NPAS4 gene expression 11 0.0047523 0.4916977 0.0180117
1028 Signaling by BMP 14 0.0021691 0.4733573 0.0095242
1123 Synthesis of PIPs at the early endosome membrane 15 0.0021666 0.4573782 0.0095242
1098 Signaling by cytosolic FGFR1 fusion mutants 17 0.0017685 0.4381012 0.0080710
1125 Synthesis of PIPs at the plasma membrane 44 0.0000013 0.4213943 0.0000096
898 Regulation of Expression and Function of Type II Classical Cadherins 16 0.0036227 0.4201913 0.0145462
901 Regulation of Homotypic Cell-Cell Adhesion 16 0.0036227 0.4201913 0.0145462
871 RND3 GTPase cycle 26 0.0003806 0.4027236 0.0019794
1199 Transcriptional Regulation by NPAS4 24 0.0007405 0.3980226 0.0036191 
head(subset(mres1$enrichment_result,s.dist<0),20) %>%
  kbl(caption="pathways downregulated in latent compared to bystander") %>%
  kable_paper("hover",full_width=FALSE)
pathways downregulated in latent compared to bystander
set setSize pANOVA s.dist p.adjustANOVA
317 Eukaryotic Translation Termination 90 0.00e+00 -0.7990742 0.0000000
1259 Viral mRNA Translation 87 0.00e+00 -0.7936818 0.0000000
315 Eukaryotic Translation Elongation 90 0.00e+00 -0.7831812 0.0000000
745 Peptide chain elongation 87 0.00e+00 -0.7795205 0.0000000
353 Formation of a pool of free 40S subunits 98 0.00e+00 -0.7728446 0.0000000
986 SRP-dependent cotranslational protein targeting to membrane 109 0.00e+00 -0.7706795 0.0000000
970 SARS-CoV-1 modulates host translation machinery 36 0.00e+00 -0.7665729 0.0000000
1006 Selenocysteine synthesis 90 0.00e+00 -0.7614834 0.0000000
344 Folding of actin by CCT/TriC 10 3.37e-05 -0.7573062 0.0002054
687 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 93 0.00e+00 -0.7461013 0.0000000
976 SARS-CoV-2 modulates host translation machinery 47 0.00e+00 -0.7457787 0.0000000
345 Formation of ATP by chemiosmotic coupling 18 1.00e-07 -0.7353620 0.0000005
394 GTP hydrolysis and joining of the 60S ribosomal subunit 109 0.00e+00 -0.7240184 0.0000000
360 Formation of the ternary complex, and subsequently, the 43S complex 50 0.00e+00 -0.7218294 0.0000000
1005 Selenoamino acid metabolism 101 0.00e+00 -0.7139292 0.0000000
541 L13a-mediated translational silencing of Ceruloplasmin expression 108 0.00e+00 -0.7137503 0.0000000
192 Cooperation of Prefoldin and TriC/CCT in actin and tubulin folding 25 0.00e+00 -0.7113725 0.0000000
955 Response of EIF2AK4 (GCN2) to amino acid deficiency 98 0.00e+00 -0.7112530 0.0000000
927 Regulation of activated PAK-2p34 by proteasome mediated degradation 48 0.00e+00 -0.7084137 0.0000000
930 Regulation of expression of SLITs and ROBOs 155 0.00e+00 -0.7072485 0.0000000 
# focus on apoptosis
mres1$enrichment_result[grep("apop",mres1$enrichment_result$set),] %>%
  kbl(caption="apoptosis pathways in latent compared to bystander") %>%
  kable_paper("hover",full_width=FALSE)
apoptosis pathways in latent compared to bystander
set setSize pANOVA s.dist p.adjustANOVA
133 Caspase activation via extrinsic apoptotic signalling pathway 21 0.1265246 0.1927006 0.2590267
201 Cytochrome c-mediated apoptotic response 12 0.3159844 0.1672310 0.4863258
1157 TP53 regulates transcription of several additional cell death genes whose specific roles in p53-dependent apoptosis remain uncertain 11 0.4754008 0.1243178 0.6280702
354 Formation of apoptosome 10 0.5635228 0.1055174 0.6996301
942 Regulation of the apoptosome activity 10 0.5635228 0.1055174 0.6996301
1145 TNFR1-induced proapoptotic signaling 23 0.6001937 0.0631701 0.7326308 
if (!file.exists("latent_v_by_mitch.html")) {
  mitch_report(res=mres1,outfile="latent_v_by_mitch.html",overwrite=TRUE)
}

# custom apop list
a1 <- m1[which(rownames(m1) %in% apop),,drop=FALSE]
a1[order(-a1$x),,drop=FALSE]
##                    x
## PTEN     21.32051053
## ST6GAL1  19.92884203
## BCL2L11   8.61209030
## MEF2C     7.76316507
## BCL2      6.86792163
## CASP9     6.83892720
## CASP3     6.44055866
## CASP7     5.63773348
## CASP8     5.56686954
## HRK       5.02172480
## SIRT1     4.48002582
## BCL2L2    4.14001595
## MCL1      3.77807893
## BCL2L1    3.49356694
## BCL2L13   1.55612052
## BBC3      0.62738833
## CDKN2A   -0.02386415
## PLEKHO2  -0.87675679
## CCR5     -1.13990223
## BAK1     -1.17457054
## BAD      -2.92209762
## BID      -3.46668487
## SELENOS  -3.71681388
## CYCS    -19.97075249
## BAX     -30.76251252
ares1 <- mitch_calc(m1,genesets=list("custom_apop"=apop))
## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.
ares1$enrichment_result
##           set setSize    pANOVA      s.dist p.adjustANOVA
## 1 custom_apop      25 0.8522829 -0.02152831     0.8522829

2. Latent (case) versus productive (ctrl).

lat_v_prod

m2 <- mitch_import(lat_v_prod,DEtype="seurat")
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 11129
## Note: no. genes in output = 11129
## Note: estimated proportion of input genes in output = 1
mres2 <- mitch_calc(x=m2,genesets=gset,cores=CORES,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
head(subset(mres2$enrichment_result,s.dist>0),20) %>%
  kbl(caption="pathways upregulated in latent as compared to productive") %>%
  kable_paper("hover",full_width=FALSE)
pathways upregulated in latent as compared to productive
set setSize pANOVA s.dist p.adjustANOVA
316 Erythropoietin activates Phosphoinositide-3-kinase (PI3K) 11 0.0002434 0.6388984 0.0019655
918 Regulation of PTEN mRNA translation 11 0.0016236 0.5488397 0.0103090
1306 tRNA processing in the mitochondrion 18 0.0000602 0.5462555 0.0005668
182 Competing endogenous RNAs (ceRNAs) regulate PTEN translation 10 0.0031086 0.5400306 0.0173024
661 NRAGE signals death through JNK 36 0.0000000 0.5385679 0.0000003
440 HDMs demethylate histones 21 0.0000260 0.5302828 0.0002611
660 NR1H3 & NR1H2 regulate gene expression linked to cholesterol transport and efflux 34 0.0000003 0.5110675 0.0000034
327 FASTK family proteins regulate processing and stability of mitochondrial RNAs 17 0.0003603 0.4998518 0.0027889
903 Regulation of CDH11 Expression and Function 15 0.0009837 0.4914762 0.0066327
915 Regulation of NPAS4 gene expression 11 0.0051956 0.4866964 0.0260375
641 Myogenesis 15 0.0011318 0.4855738 0.0074017
1106 Signaling by cytosolic FGFR1 fusion mutants 17 0.0006134 0.4799792 0.0044820
905 Regulation of Expression and Function of Type II Classical Cadherins 16 0.0010964 0.4714748 0.0072430
908 Regulation of Homotypic Cell-Cell Adhesion 16 0.0010964 0.4714748 0.0072430
312 Endosomal/Vacuolar pathway 10 0.0107292 0.4660131 0.0457125
659 NR1H2 and NR1H3-mediated signaling 39 0.0000007 0.4595593 0.0000085
1232 Translocation of ZAP-70 to Immunological synapse 10 0.0131426 0.4529544 0.0527315
334 FGFR1 mutant receptor activation 20 0.0005644 0.4454856 0.0041711
526 Interleukin-6 signaling 10 0.0149014 0.4447163 0.0575054
9 ALK mutants bind TKIs 12 0.0083224 0.4400018 0.0365291 
head(subset(mres2$enrichment_result,s.dist<0),20) %>%
  kbl(caption="pathways downregulated in latent as compared to productive") %>%
  kable_paper("hover",full_width=FALSE)
pathways downregulated in latent as compared to productive
set setSize pANOVA s.dist p.adjustANOVA
18 APOBEC3G mediated resistance to HIV-1 infection 10 0.0000265 -0.7671913 0.0002629
624 Mitochondrial translation termination 88 0.0000000 -0.6680328 0.0000000
622 Mitochondrial translation elongation 88 0.0000000 -0.6673803 0.0000000
623 Mitochondrial translation initiation 88 0.0000000 -0.6654103 0.0000000
621 Mitochondrial translation 94 0.0000000 -0.6547369 0.0000000
599 Metabolism of polyamines 54 0.0000000 -0.6408227 0.0000000
946 Regulation of ornithine decarboxylase (ODC) 49 0.0000000 -0.6285641 0.0000000
98 Binding and entry of HIV virion 11 0.0004553 -0.6104761 0.0034625
349 Folding of actin by CCT/TriC 10 0.0009415 -0.6040471 0.0063805
1264 Vif-mediated degradation of APOBEC3G 53 0.0000000 -0.5915902 0.0000000
1115 Somitogenesis 46 0.0000000 -0.5856364 0.0000000
328 FBXL7 down-regulates AURKA during mitotic entry and in early mitosis 52 0.0000000 -0.5849300 0.0000000
675 Negative regulation of NOTCH4 signaling 53 0.0000000 -0.5769844 0.0000000
934 Regulation of activated PAK-2p34 by proteasome mediated degradation 48 0.0000000 -0.5740344 0.0000000
619 Mitochondrial protein import 63 0.0000000 -0.5725244 0.0000000
199 Cross-presentation of soluble exogenous antigens (endosomes) 46 0.0000000 -0.5701878 0.0000000
1254 Ubiquitin-dependent degradation of Cyclin D 50 0.0000000 -0.5636971 0.0000000
1253 Ubiquitin Mediated Degradation of Phosphorylated Cdc25A 49 0.0000000 -0.5622596 0.0000000
1293 p53-Independent DNA Damage Response 49 0.0000000 -0.5622596 0.0000000
1294 p53-Independent G1/S DNA damage checkpoint 49 0.0000000 -0.5622596 0.0000000 
# focus on apoptosis
mres2$enrichment_result[grep("apop",mres2$enrichment_result$set),] %>%
  kbl(caption="apoptosis pathways in latent compared to productive") %>%
  kable_paper("hover",full_width=FALSE)
apoptosis pathways in latent compared to productive
set setSize pANOVA s.dist p.adjustANOVA
136 Caspase activation via extrinsic apoptotic signalling pathway 21 0.1559307 0.1789444 0.3209608
1165 TP53 regulates transcription of several additional cell death genes whose specific roles in p53-dependent apoptosis remain uncertain 11 0.5955956 0.0924463 0.7645133
1153 TNFR1-induced proapoptotic signaling 23 0.6692313 0.0514880 0.8108062
205 Cytochrome c-mediated apoptotic response 12 0.7745703 0.0477647 0.8809895
359 Formation of apoptosome 10 0.9685042 -0.0072129 0.9807009
949 Regulation of the apoptosome activity 10 0.9685042 -0.0072129 0.9807009 
if (!file.exists("latent_v_prod_mitch.html")) {
  mitch_report(res=mres2,outfile="latent_v_prod_mitch.html",overwrite=TRUE)
}

# custom apop list
a2 <- m2[which(rownames(m2) %in% apop),,drop=FALSE]
a2[order(-a2$x),,drop=FALSE]
##                    x
## ST6GAL1  13.53052980
## BCL2L11  10.56552806
## PTEN      4.92530925
## MCL1      4.27711111
## CASP9     4.15765685
## CASP3     4.08280387
## HRK       3.18041204
## PLEKHO2   1.33181829
## SIRT1     1.20325251
## CASP8     1.13464266
## CDKN2A    1.11444415
## CCR5      0.85072500
## CASP7     0.68857954
## BCL2L1    0.25803111
## BCL2L2    0.25370945
## MEF2C     0.19693914
## BCL2      0.02468332
## BBC3     -0.39134438
## BCL2L13  -1.53199007
## BAK1     -1.67359684
## SELENOS  -3.45887483
## BID      -9.77445785
## BAD     -10.94394139
## BAX     -11.92870938
## CYCS    -31.55237715
ares2 <- mitch_calc(m2,genesets=list("custom_apop"=apop))
## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.
ares2$enrichment_result
##           set setSize   pANOVA      s.dist p.adjustANOVA
## 1 custom_apop      25 0.691536 -0.04587176      0.691536

3. Productive (case) vs all other clusters (ctrl).

prod_v_all

m3 <- mitch_import(prod_v_all,DEtype="seurat")
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 10703
## Note: no. genes in output = 10703
## Note: estimated proportion of input genes in output = 1
mres3 <- mitch_calc(x=m3,genesets=gset,cores=CORES,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
head(subset(mres3$enrichment_result,s.dist>0),20) %>%
  kbl(caption="pathways upregulated in productive as compared to all other cell clusters") %>%
  kable_paper("hover",full_width=FALSE)
pathways upregulated in productive as compared to all other cell clusters
set setSize pANOVA s.dist p.adjustANOVA
158 Chemokine receptors bind chemokines 10 0.0016100 0.5760591 0.0175472
1275 mRNA decay by 3’ to 5’ exoribonuclease 13 0.0026381 0.4817587 0.0259210
18 APOBEC3G mediated resistance to HIV-1 infection 10 0.0100276 0.4703077 0.0710700
589 Metabolism of folate and pterines 14 0.0041997 0.4419898 0.0373081
1 2-LTR circle formation 12 0.0120528 0.4186543 0.0811532
785 Processing and activation of SUMO 10 0.0248828 0.4097447 0.1407099
97 Binding and entry of HIV virion 11 0.0194416 0.4070163 0.1156685
10 APC truncation mutants have impaired AXIN binding 12 0.0369386 0.3479406 0.1751986
24 AXIN missense mutants destabilize the destruction complex 12 0.0369386 0.3479406 0.1751986
1026 Signaling by AMER1 mutants 12 0.0369386 0.3479406 0.1751986
1027 Signaling by APC mutants 12 0.0369386 0.3479406 0.1751986
1028 Signaling by AXIN mutants 12 0.0369386 0.3479406 0.1751986
1240 Truncations of AMER1 destabilize the destruction complex 12 0.0369386 0.3479406 0.1751986
500 Integration of provirus 14 0.0354679 0.3247130 0.1751986
505 Interconversion of nucleotide di- and triphosphates 24 0.0072996 0.3165168 0.0564554
481 Inactivation, recovery and regulation of the phototransduction cascade 16 0.0299433 0.3135702 0.1585163
256 Disassembly of the destruction complex and recruitment of AXIN to the membrane 18 0.0221585 0.3115687 0.1278173
1248 VEGFR2 mediated cell proliferation 15 0.0386276 0.3085205 0.1795697
1174 The phototransduction cascade 17 0.0368939 0.2924662 0.1751986
1279 p38MAPK events 11 0.0967718 0.2892562 0.3125264 
head(subset(mres3$enrichment_result,s.dist<0),20) %>%
  kbl(caption="pathways downregulated in productive as compared to all other cell clusters") %>%
  kable_paper("hover",full_width=FALSE)
pathways downregulated in productive as compared to all other cell clusters
set setSize pANOVA s.dist p.adjustANOVA
318 Eukaryotic Translation Elongation 90 0.00e+00 -0.7595119 0.0000000
746 Peptide chain elongation 87 0.00e+00 -0.7570908 0.0000000
1256 Viral mRNA Translation 87 0.00e+00 -0.7430846 0.0000000
355 Formation of a pool of free 40S subunits 98 0.00e+00 -0.7329956 0.0000000
1222 Translocation of ZAP-70 to Immunological synapse 10 6.39e-05 -0.7300664 0.0011201
971 SARS-CoV-1 modulates host translation machinery 36 0.00e+00 -0.7203551 0.0000000
542 L13a-mediated translational silencing of Ceruloplasmin expression 108 0.00e+00 -0.6969229 0.0000000
396 GTP hydrolysis and joining of the 60S ribosomal subunit 109 0.00e+00 -0.6936807 0.0000000
1007 Selenocysteine synthesis 90 0.00e+00 -0.6841222 0.0000000
987 SRP-dependent cotranslational protein targeting to membrane 109 0.00e+00 -0.6726440 0.0000000
688 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 93 0.00e+00 -0.6719011 0.0000000
398 Gap junction assembly 10 2.99e-04 -0.6603947 0.0042159
1230 Transport of connexons to the plasma membrane 10 2.99e-04 -0.6603947 0.0042159
320 Eukaryotic Translation Termination 91 0.00e+00 -0.6590766 0.0000000
956 Response of EIF2AK4 (GCN2) to amino acid deficiency 98 0.00e+00 -0.6559873 0.0000000
408 Generation of second messenger molecules 18 1.60e-06 -0.6540633 0.0000341
127 Cap-dependent Translation Initiation 116 0.00e+00 -0.6480386 0.0000000
319 Eukaryotic Translation Initiation 116 0.00e+00 -0.6480386 0.0000000
362 Formation of the ternary complex, and subsequently, the 43S complex 50 0.00e+00 -0.6359467 0.0000000
756 Phosphorylation of CD3 and TCR zeta chains 14 5.50e-05 -0.6225759 0.0009913 
# focus on apoptosis
mres3$enrichment_result[grep("apop",mres3$enrichment_result$set),] %>%
  kbl(caption="apoptosis pathways in productive compared to all other cell clusters") %>%
  kable_paper("hover",full_width=FALSE)
apoptosis pathways in productive compared to all other cell clusters
set setSize pANOVA s.dist p.adjustANOVA
356 Formation of apoptosome 10 0.2732745 0.2001309 0.5389376
943 Regulation of the apoptosome activity 10 0.2732745 0.2001309 0.5389376
204 Cytochrome c-mediated apoptotic response 12 0.2636585 0.1864185 0.5297198
1146 TNFR1-induced proapoptotic signaling 23 0.3913076 0.1033219 0.6432522
1156 TP53 regulates transcription of several additional cell death genes whose specific roles in p53-dependent apoptosis remain uncertain 10 0.6132225 0.0923408 0.7961457
135 Caspase activation via extrinsic apoptotic signalling pathway 20 0.5637497 0.0746045 0.7652250 
if (!file.exists("prod_v_all_mitch.html")) {
  mitch_report(res=mres3,outfile="prod_v_all_mitch.html",overwrite=TRUE)
}

# custom apop list
a3 <- m3[which(rownames(m3) %in% apop),,drop=FALSE]
a3[order(-a3$x),,drop=FALSE]
##                    x
## PTEN     14.92296985
## BCL2     14.50872197
## MEF2C    14.10634431
## BCL2L13  12.43608061
## BAD      11.99443548
## BID       9.38496521
## CYCS      9.36009013
## CASP7     6.22697375
## BCL2L1    5.49833946
## BCL2L2    5.28215310
## BBC3      3.96761906
## CASP8     3.77896500
## SELENOS  -0.03685491
## CCR5     -0.03688655
## HRK      -0.12479074
## CASP9    -0.26071214
## CASP3    -0.41067524
## ST6GAL1  -0.68162839
## BAK1     -1.09534842
## MCL1     -1.47518417
## SIRT1    -1.78976622
## CDKN2A   -2.58024064
## BCL2L11  -3.29009416
## BAX      -9.03875061
## PLEKHO2 -11.74504340
ares3 <- mitch_calc(m3,genesets=list("custom_apop"=apop))
## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.
ares3$enrichment_result
##           set setSize    pANOVA     s.dist p.adjustANOVA
## 1 custom_apop      25 0.6553037 0.05161266     0.6553037

4. Productive (case) vs mock (ctrl)

prod_v_mock

m4 <- mitch_import(prod_v_mock,DEtype="seurat")
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 10713
## Note: no. genes in output = 10713
## Note: estimated proportion of input genes in output = 1
mres4 <- mitch_calc(x=m4,genesets=gset,cores=CORES,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
head(subset(mres4$enrichment_result,s.dist>0),20) %>%
  kbl(caption="pathways upregulated in productive as compared to mock-infected") %>%
  kable_paper("hover",full_width=FALSE)
pathways upregulated in productive as compared to mock-infected
set setSize pANOVA s.dist p.adjustANOVA
158 Chemokine receptors bind chemokines 10 0.0006616 0.6218630 0.0080975
1273 mRNA decay by 3’ to 5’ exoribonuclease 13 0.0005141 0.5564198 0.0064636
18 APOBEC3G mediated resistance to HIV-1 infection 10 0.0046658 0.5167523 0.0351290
97 Binding and entry of HIV virion 11 0.0179723 0.4121065 0.0908421
1 2-LTR circle formation 12 0.0295174 0.3629567 0.1239547
588 Metabolism of folate and pterines 14 0.0299192 0.3352383 0.1249852
612 Mitochondrial iron-sulfur cluster biogenesis 13 0.0404781 0.3282818 0.1583659
784 Processing and activation of SUMO 10 0.0846328 0.3149958 0.2469055
258 Diseases associated with N-glycosylation of proteins 20 0.0150291 0.3142149 0.0787963
593 Metabolism of polyamines 54 0.0000681 0.3136282 0.0012424
134 Caspase activation via Death Receptors in the presence of ligand 14 0.0451757 0.3092812 0.1690823
536 KSRP (KHSRP) binds and destabilizes mRNA 16 0.0354422 0.3038118 0.1434303
1237 Tristetraprolin (TTP, ZFP36) binds and destabilizes mRNA 17 0.0329030 0.2989573 0.1348398
939 Regulation of ornithine decarboxylase (ODC) 49 0.0004208 0.2914938 0.0053952
884 Receptor Mediated Mitophagy 11 0.0951539 0.2906678 0.2661431
1150 TP53 Regulates Transcription of Cell Death Genes 36 0.0038017 0.2789381 0.0306345
586 Metabolism of cofactors 16 0.0577821 0.2740722 0.1958993
1246 VEGFR2 mediated cell proliferation 15 0.0694315 0.2708544 0.2209184
361 Formation of the cornified envelope 14 0.0905623 0.2613328 0.2603593
539 Keratinization 14 0.0905623 0.2613328 0.2603593 
head(subset(mres4$enrichment_result,s.dist<0),20) %>%
  kbl(caption="pathways downregulated in productive as compared to mock-infected") %>%
  kable_paper("hover",full_width=FALSE)
pathways downregulated in productive as compared to mock-infected
set setSize pANOVA s.dist p.adjustANOVA
745 Peptide chain elongation 87 0.0000000 -0.7112180 0.0000000
318 Eukaryotic Translation Elongation 90 0.0000000 -0.7043982 0.0000000
1254 Viral mRNA Translation 87 0.0000000 -0.6986377 0.0000000
1221 Translocation of ZAP-70 to Immunological synapse 10 0.0001366 -0.6965524 0.0021059
355 Formation of a pool of free 40S subunits 98 0.0000000 -0.6847914 0.0000000
970 SARS-CoV-1 modulates host translation machinery 36 0.0000000 -0.6779058 0.0000000
541 L13a-mediated translational silencing of Ceruloplasmin expression 108 0.0000000 -0.6525905 0.0000000
396 GTP hydrolysis and joining of the 60S ribosomal subunit 109 0.0000000 -0.6457006 0.0000000
1006 Selenocysteine synthesis 90 0.0000000 -0.6394908 0.0000000
407 Generation of second messenger molecules 18 0.0000027 -0.6388343 0.0000673
325 FASTK family proteins regulate processing and stability of mitochondrial RNAs 17 0.0000052 -0.6384245 0.0001266
687 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 93 0.0000000 -0.6296357 0.0000000
986 SRP-dependent cotranslational protein targeting to membrane 109 0.0000000 -0.6178281 0.0000000
320 Eukaryotic Translation Termination 91 0.0000000 -0.6175034 0.0000000
1293 tRNA processing in the mitochondrion 18 0.0000069 -0.6121552 0.0001659
755 Phosphorylation of CD3 and TCR zeta chains 14 0.0000742 -0.6116593 0.0013163
955 Response of EIF2AK4 (GCN2) to amino acid deficiency 98 0.0000000 -0.6102012 0.0000000
127 Cap-dependent Translation Initiation 116 0.0000000 -0.6072652 0.0000000
319 Eukaryotic Translation Initiation 116 0.0000000 -0.6072652 0.0000000
362 Formation of the ternary complex, and subsequently, the 43S complex 50 0.0000000 -0.5892826 0.0000000 
# focus on apoptosis
mres4$enrichment_result[grep("apop",mres4$enrichment_result$set),] %>%
  kbl(caption="apoptosis pathways in productive compared to mock-infected") %>%
  kable_paper("hover",full_width=FALSE)
apoptosis pathways in productive compared to mock-infected
set setSize pANOVA s.dist p.adjustANOVA
356 Formation of apoptosome 10 0.2139934 0.2270018 0.4384833
942 Regulation of the apoptosome activity 10 0.2139934 0.2270018 0.4384833
204 Cytochrome c-mediated apoptotic response 12 0.3215775 0.1653117 0.5394338
1155 TP53 regulates transcription of several additional cell death genes whose specific roles in p53-dependent apoptosis remain uncertain 10 0.6824221 0.0747454 0.8337138
135 Caspase activation via extrinsic apoptotic signalling pathway 20 0.6112032 0.0656972 0.7881025
1145 TNFR1-induced proapoptotic signaling 23 0.7249732 0.0424045 0.8566061 
if (!file.exists("prod_v_mock_mitch.html")) {
  mitch_report(res=mres4,outfile="prod_v_mock_mitch.html",overwrite=TRUE)
}

# custom apop list
a4 <- m4[which(rownames(m4) %in% apop),,drop=FALSE]
a4[order(-a4$x),,drop=FALSE]
##                    x
## PTEN     11.27138605
## BAD       8.90529930
## BID       8.71811788
## BCL2      8.12345445
## BCL2L13   7.85896810
## MEF2C     6.56896458
## CYCS      4.74623047
## CASP7     4.66258312
## BBC3      4.09279560
## BCL2L1    3.58884604
## BCL2L2    2.47850099
## BAK1      1.59402864
## HRK      -0.09151097
## CCR5     -0.17021668
## CASP9    -0.33859942
## SELENOS  -0.46553562
## SIRT1    -0.64886280
## CDKN2A   -1.32598799
## ST6GAL1  -1.37049203
## CASP8    -1.85377861
## MCL1     -2.22019479
## BCL2L11  -3.82304227
## CASP3    -4.04225029
## BAX      -5.32460664
## PLEKHO2 -16.28613144
ares4 <- mitch_calc(m4,genesets=list("custom_apop"=apop))
## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.
ares4$enrichment_result
##           set setSize   pANOVA      s.dist p.adjustANOVA
## 1 custom_apop      25 0.917809 -0.01193114      0.917809

5. Latent (case) vs mock (ctrl)

lat_v_mock

m5 <- mitch_import(prod_v_mock,DEtype="seurat")
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 10713
## Note: no. genes in output = 10713
## Note: estimated proportion of input genes in output = 1
mres5 <- mitch_calc(x=m5,genesets=gset,cores=CORES,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
head(subset(mres5$enrichment_result,s.dist>0),20) %>%
  kbl(caption="pathways upregulated in latent as compared to mock-infected") %>%
  kable_paper("hover",full_width=FALSE)
pathways upregulated in latent as compared to mock-infected
set setSize pANOVA s.dist p.adjustANOVA
158 Chemokine receptors bind chemokines 10 0.0006616 0.6218630 0.0080975
1273 mRNA decay by 3’ to 5’ exoribonuclease 13 0.0005141 0.5564198 0.0064636
18 APOBEC3G mediated resistance to HIV-1 infection 10 0.0046658 0.5167523 0.0351290
97 Binding and entry of HIV virion 11 0.0179723 0.4121065 0.0908421
1 2-LTR circle formation 12 0.0295174 0.3629567 0.1239547
588 Metabolism of folate and pterines 14 0.0299192 0.3352383 0.1249852
612 Mitochondrial iron-sulfur cluster biogenesis 13 0.0404781 0.3282818 0.1583659
784 Processing and activation of SUMO 10 0.0846328 0.3149958 0.2469055
258 Diseases associated with N-glycosylation of proteins 20 0.0150291 0.3142149 0.0787963
593 Metabolism of polyamines 54 0.0000681 0.3136282 0.0012424
134 Caspase activation via Death Receptors in the presence of ligand 14 0.0451757 0.3092812 0.1690823
536 KSRP (KHSRP) binds and destabilizes mRNA 16 0.0354422 0.3038118 0.1434303
1237 Tristetraprolin (TTP, ZFP36) binds and destabilizes mRNA 17 0.0329030 0.2989573 0.1348398
939 Regulation of ornithine decarboxylase (ODC) 49 0.0004208 0.2914938 0.0053952
884 Receptor Mediated Mitophagy 11 0.0951539 0.2906678 0.2661431
1150 TP53 Regulates Transcription of Cell Death Genes 36 0.0038017 0.2789381 0.0306345
586 Metabolism of cofactors 16 0.0577821 0.2740722 0.1958993
1246 VEGFR2 mediated cell proliferation 15 0.0694315 0.2708544 0.2209184
361 Formation of the cornified envelope 14 0.0905623 0.2613328 0.2603593
539 Keratinization 14 0.0905623 0.2613328 0.2603593 
head(subset(mres5$enrichment_result,s.dist<0),20) %>%
  kbl(caption="pathways downregulated in latent as compared to mock-infected") %>%
  kable_paper("hover",full_width=FALSE)
pathways downregulated in latent as compared to mock-infected
set setSize pANOVA s.dist p.adjustANOVA
745 Peptide chain elongation 87 0.0000000 -0.7112180 0.0000000
318 Eukaryotic Translation Elongation 90 0.0000000 -0.7043982 0.0000000
1254 Viral mRNA Translation 87 0.0000000 -0.6986377 0.0000000
1221 Translocation of ZAP-70 to Immunological synapse 10 0.0001366 -0.6965524 0.0021059
355 Formation of a pool of free 40S subunits 98 0.0000000 -0.6847914 0.0000000
970 SARS-CoV-1 modulates host translation machinery 36 0.0000000 -0.6779058 0.0000000
541 L13a-mediated translational silencing of Ceruloplasmin expression 108 0.0000000 -0.6525905 0.0000000
396 GTP hydrolysis and joining of the 60S ribosomal subunit 109 0.0000000 -0.6457006 0.0000000
1006 Selenocysteine synthesis 90 0.0000000 -0.6394908 0.0000000
407 Generation of second messenger molecules 18 0.0000027 -0.6388343 0.0000673
325 FASTK family proteins regulate processing and stability of mitochondrial RNAs 17 0.0000052 -0.6384245 0.0001266
687 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 93 0.0000000 -0.6296357 0.0000000
986 SRP-dependent cotranslational protein targeting to membrane 109 0.0000000 -0.6178281 0.0000000
320 Eukaryotic Translation Termination 91 0.0000000 -0.6175034 0.0000000
1293 tRNA processing in the mitochondrion 18 0.0000069 -0.6121552 0.0001659
755 Phosphorylation of CD3 and TCR zeta chains 14 0.0000742 -0.6116593 0.0013163
955 Response of EIF2AK4 (GCN2) to amino acid deficiency 98 0.0000000 -0.6102012 0.0000000
127 Cap-dependent Translation Initiation 116 0.0000000 -0.6072652 0.0000000
319 Eukaryotic Translation Initiation 116 0.0000000 -0.6072652 0.0000000
362 Formation of the ternary complex, and subsequently, the 43S complex 50 0.0000000 -0.5892826 0.0000000 
# focus on apoptosis
mres5$enrichment_result[grep("apop",mres5$enrichment_result$set),] %>%
  kbl(caption="apoptosis pathways in latent compared to mock-infected") %>%
  kable_paper("hover",full_width=FALSE)
apoptosis pathways in latent compared to mock-infected
set setSize pANOVA s.dist p.adjustANOVA
356 Formation of apoptosome 10 0.2139934 0.2270018 0.4384833
942 Regulation of the apoptosome activity 10 0.2139934 0.2270018 0.4384833
204 Cytochrome c-mediated apoptotic response 12 0.3215775 0.1653117 0.5394338
1155 TP53 regulates transcription of several additional cell death genes whose specific roles in p53-dependent apoptosis remain uncertain 10 0.6824221 0.0747454 0.8337138
135 Caspase activation via extrinsic apoptotic signalling pathway 20 0.6112032 0.0656972 0.7881025
1145 TNFR1-induced proapoptotic signaling 23 0.7249732 0.0424045 0.8566061 
if (!file.exists("latent_v_mock_mitch.html")) {
  mitch_report(res=mres5,outfile="latent_v_mock_mitch.html",overwrite=TRUE)
}

# custom apop list
a5 <- m5[which(rownames(m5) %in% apop),,drop=FALSE]
a5[order(-a5$x),,drop=FALSE]
##                    x
## PTEN     11.27138605
## BAD       8.90529930
## BID       8.71811788
## BCL2      8.12345445
## BCL2L13   7.85896810
## MEF2C     6.56896458
## CYCS      4.74623047
## CASP7     4.66258312
## BBC3      4.09279560
## BCL2L1    3.58884604
## BCL2L2    2.47850099
## BAK1      1.59402864
## HRK      -0.09151097
## CCR5     -0.17021668
## CASP9    -0.33859942
## SELENOS  -0.46553562
## SIRT1    -0.64886280
## CDKN2A   -1.32598799
## ST6GAL1  -1.37049203
## CASP8    -1.85377861
## MCL1     -2.22019479
## BCL2L11  -3.82304227
## CASP3    -4.04225029
## BAX      -5.32460664
## PLEKHO2 -16.28613144
ares5 <- mitch_calc(m5,genesets=list("custom_apop"=apop))
## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.
ares5$enrichment_result
##           set setSize   pANOVA      s.dist p.adjustANOVA
## 1 custom_apop      25 0.917809 -0.01193114      0.917809

6. Multidimensional

ml <- list("lat_v_by"= lat_v_by,
  "lat_v_prod" = lat_v_prod,
  "prod_v_all" = prod_v_all,
  "prod_v_mock" = prod_v_mock,
  "lat_v_mock" = lat_v_mock)

str(ml)
## List of 5
##  $ lat_v_by   :'data.frame': 10939 obs. of  5 variables:
##   ..$ p_val     : num [1:10939] 0.00 1.07e-175 5.21e-162 1.99e-143 5.99e-125 ...
##   ..$ avg_log2FC: num [1:10939] 0.777 0.371 0.272 0.316 0.353 ...
##   ..$ pct.1     : num [1:10939] 0.987 0.599 0.434 0.484 0.582 0.198 1 0.998 0.133 1 ...
##   ..$ pct.2     : num [1:10939] 0.327 0.121 0.06 0.09 0.155 0.018 1 1 0.008 0.999 ...
##   ..$ p_val_adj : num [1:10939] 0.00 3.93e-171 1.91e-157 7.29e-139 2.19e-120 ...
##  $ lat_v_prod :'data.frame': 11129 obs. of  5 variables:
##   ..$ p_val     : num [1:11129] 3.27e-213 2.00e-197 1.11e-195 1.30e-182 2.45e-182 ...
##   ..$ avg_log2FC: num [1:11129] -1.37 -1.07 -1.25 -1.09 -1.03 ...
##   ..$ pct.1     : num [1:11129] 0.484 0.582 0.599 0.987 0.434 0.198 0.133 0.061 1 1 ...
##   ..$ pct.2     : num [1:11129] 0.956 0.953 0.946 0.971 0.922 0.8 0.741 0.647 1 1 ...
##   ..$ p_val_adj : num [1:11129] 1.20e-208 7.34e-193 4.06e-191 4.75e-178 8.96e-178 ...
##  $ prod_v_all :'data.frame': 10703 obs. of  5 variables:
##   ..$ p_val     : num [1:10703] 0 0 0 0 0 0 0 0 0 0 ...
##   ..$ avg_log2FC: num [1:10703] 0.0807 -0.1934 -0.2092 1.5865 1.6545 ...
##   ..$ pct.1     : num [1:10703] 1 1 1 0.946 0.956 0.971 0.8 0.498 0.922 0.647 ...
##   ..$ pct.2     : num [1:10703] 1 1 1 0.156 0.124 0.366 0.033 0.009 0.086 0.016 ...
##   ..$ p_val_adj : num [1:10703] 0 0 0 0 0 0 0 0 0 0 ...
##  $ prod_v_mock:'data.frame': 10713 obs. of  5 variables:
##   ..$ p_val     : num [1:10713] 0 0 0 0 0 ...
##   ..$ avg_log2FC: num [1:10713] 1.627 1.684 1.886 0.806 1.314 ...
##   ..$ pct.1     : num [1:10713] 0.946 0.956 0.971 0.8 0.922 0.647 0.741 0.953 0.603 1 ...
##   ..$ pct.2     : num [1:10713] 0.101 0.092 0.269 0.019 0.039 0.015 0.009 0.082 0.004 0.999 ...
##   ..$ p_val_adj : num [1:10713] 0 0 0 0 0 ...
##  $ lat_v_mock :'data.frame': 10951 obs. of  5 variables:
##   ..$ p_val     : num [1:10951] 1.99e-272 1.75e-152 4.60e-143 3.08e-132 2.60e-100 ...
##   ..$ avg_log2FC: num [1:10951] 0.798 0.391 0.378 0.282 0.316 ...
##   ..$ pct.1     : num [1:10951] 0.987 0.582 0.599 0.434 0.484 1 0.993 1 0.198 1 ...
##   ..$ pct.2     : num [1:10951] 0.269 0.082 0.101 0.039 0.092 0.998 0.982 1 0.019 1 ...
##   ..$ p_val_adj : num [1:10951] 7.29e-268 6.40e-148 1.68e-138 1.13e-127 9.54e-96 ...
mm <- mitch_import(ml,DEtype="seurat")
## Note: Mean no. genes in input = 10887
## Note: no. genes in output = 11135
## Note: estimated proportion of input genes in output = 1.02
str(mm)
## 'data.frame':    11135 obs. of  5 variables:
##  $ lat_v_by   : num  175 175 161 143 124 ...
##  $ lat_v_prod : num  -182 -195 -182 -212 -197 ...
##  $ prod_v_all : num  221 221 221 221 221 ...
##  $ prod_v_mock: num  261 261 261 261 261 ...
##  $ lat_v_mock : num  271.7 142.3 131.5 99.6 151.8 ...
head(mm)
##      lat_v_by lat_v_prod prod_v_all prod_v_mock lat_v_mock
## pol 174.97945  -181.8869    221.219    260.9277  271.70110
## env 174.96945  -194.9551    221.219    260.9277  142.33738
## tat 161.28311  -181.6113    221.219    260.9277  131.51163
## gag 142.70074  -212.4848    221.219    260.9277   99.58430
## nef 124.22267  -196.6982    221.219    260.9277  151.75740
## vif  90.52759  -139.3709    221.219    260.9277   53.39414
tail(mm)
##            lat_v_by  lat_v_prod prod_v_all prod_v_mock lat_v_mock
## SMC5-AS1         NA  0.11864561  10.942156    6.308547         NA
## AL157938.2       NA -0.08461744   5.984703    4.721785         NA
## BX537318.1       NA  0.07932208   2.314550    3.143039         NA
## PCLAF            NA          NA -29.005164  -40.016357 -3.5194260
## LINC01678        NA          NA         NA   -4.533175 -0.7921354
## COL22A1          NA          NA         NA   -1.824147 -4.5226730
mresm <- mitch_calc(x=mm,genesets=gset,cores=CORES,priority="effect")
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
head(mresm$enrichment_result,50) %>%
  kbl(caption="pathways dysregulated in all") %>%
  kable_paper("hover",full_width=FALSE)
pathways dysregulated in all
set setSize pMANOVA s.lat_v_by s.lat_v_prod s.prod_v_all s.prod_v_mock s.lat_v_mock p.lat_v_by p.lat_v_prod p.prod_v_all p.prod_v_mock p.lat_v_mock s.dist SD p.adjustMANOVA
1267 Viral mRNA Translation 87 0.0000000 -0.7797113 0.1712763 -0.7142555 -0.6721604 -0.6629507 0.0000000 0.0043935 0.0000000 0.0000000 0.0000000 1.427848 0.3955906 0.0000000
751 Peptide chain elongation 87 0.0000000 -0.7657992 0.1890858 -0.7277183 -0.6842639 -0.6415616 0.0000000 0.0017222 0.0000000 0.0000000 0.0000000 1.425340 0.4024747 0.0000000
320 Eukaryotic Translation Elongation 90 0.0000000 -0.7693956 0.1681593 -0.7300454 -0.6777026 -0.6447621 0.0000000 0.0044285 0.0000000 0.0000000 0.0000000 1.424153 0.3936173 0.0000000
358 Formation of a pool of free 40S subunits 98 0.0000000 -0.7592409 0.1747784 -0.7045579 -0.6588389 -0.6360596 0.0000000 0.0020592 0.0000000 0.0000000 0.0000000 1.393570 0.3894577 0.0000000
978 SARS-CoV-1 modulates host translation machinery 36 0.0000000 -0.7530796 0.0836163 -0.6924078 -0.6522142 -0.6202899 0.0000000 0.3529491 0.0000000 0.0000000 0.0000000 1.365171 0.3448570 0.0000000
322 Eukaryotic Translation Termination 91 0.0000000 -0.7850087 0.1093175 -0.6335067 -0.5941009 -0.6598258 0.0000000 0.0416377 0.0000000 0.0000000 0.0000000 1.348275 0.3549381 0.0000000
1014 Selenocysteine synthesis 90 0.0000000 -0.7480796 0.1504679 -0.6575806 -0.6152551 -0.6242130 0.0000000 0.0105557 0.0000000 0.0000000 0.0000000 1.335237 0.3668074 0.0000000
994 SRP-dependent cotranslational protein targeting to membrane 109 0.0000000 -0.7571139 0.0818686 -0.6465477 -0.5944134 -0.6394890 0.0000000 0.1139386 0.0000000 0.0000000 0.0000000 1.326733 0.3368626 0.0000000
693 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 93 0.0000000 -0.7329683 0.1423701 -0.6458337 -0.6057734 -0.6166383 0.0000000 0.0137135 0.0000000 0.0000000 0.0000000 1.312180 0.3579926 0.0000000
399 GTP hydrolysis and joining of the 60S ribosomal subunit 109 0.0000000 -0.7112741 0.1524794 -0.6667682 -0.6212295 -0.5987246 0.0000000 0.0043968 0.0000000 0.0000000 0.0000000 1.310778 0.3612573 0.0000000
546 L13a-mediated translational silencing of Ceruloplasmin expression 108 0.0000000 -0.7011868 0.1617369 -0.6698847 -0.6278583 -0.5890318 0.0000000 0.0026940 0.0000000 0.0000000 0.0000000 1.306802 0.3641492 0.0000000
350 Formation of ATP by chemiosmotic coupling 18 0.0000000 -0.7224181 -0.4827246 -0.4922313 -0.4008032 -0.7206277 0.0000001 0.0004842 0.0002121 0.0026019 0.0000001 1.295048 0.1483390 0.0000000
963 Response of EIF2AK4 (GCN2) to amino acid deficiency 98 0.0000000 -0.6987334 0.1605868 -0.6305372 -0.5870755 -0.5860275 0.0000000 0.0045136 0.0000000 0.0000000 0.0000000 1.264784 0.3545724 0.0000000
128 Cap-dependent Translation Initiation 116 0.0000000 -0.6868869 0.1219781 -0.6228969 -0.5842508 -0.5770241 0.0000000 0.0176566 0.0000000 0.0000000 0.0000000 1.244588 0.3336775 0.0000000
321 Eukaryotic Translation Initiation 116 0.0000000 -0.6868869 0.1219781 -0.6228969 -0.5842508 -0.5770241 0.0000000 0.0176566 0.0000000 0.0000000 0.0000000 1.244588 0.3336775 0.0000000
365 Formation of the ternary complex, and subsequently, the 43S complex 50 0.0000000 -0.7091236 0.0754102 -0.6112741 -0.5669496 -0.5846938 0.0000000 0.3201406 0.0000000 0.0000000 0.0000000 1.243181 0.3149324 0.0000000
1013 Selenoamino acid metabolism 101 0.0000000 -0.7013625 0.1059072 -0.5803125 -0.5419502 -0.5823026 0.0000000 0.0526854 0.0000000 0.0000000 0.0000000 1.213538 0.3219688 0.0000000
1232 Translocation of ZAP-70 to Immunological synapse 10 0.0000007 -0.2897665 0.4527103 -0.7017423 -0.6701541 -0.4457775 0.1246647 0.0122509 0.0000752 0.0001561 0.0158844 1.195482 0.4694548 0.0000049
349 Folding of actin by CCT/TriC 10 0.0000733 -0.7439760 -0.6037216 -0.2047523 -0.0520319 -0.6817182 0.0000281 0.0003503 0.2883756 0.8072446 0.0000707 1.194718 0.3090342 0.0004131
622 Mitochondrial translation elongation 88 0.0000000 -0.6921542 -0.6670206 0.0859292 0.2193056 -0.6412883 0.0000000 0.0000000 0.1093302 0.0001374 0.0000000 1.179289 0.4516535 0.0000000
624 Mitochondrial translation termination 88 0.0000000 -0.6746850 -0.6676728 0.1097298 0.2324849 -0.6256711 0.0000000 0.0000000 0.0451621 0.0000546 0.0000000 1.165565 0.4554906 0.0000000
623 Mitochondrial translation initiation 88 0.0000000 -0.6769320 -0.6650518 0.0987463 0.2219028 -0.6310600 0.0000000 0.0000000 0.0691170 0.0001150 0.0000000 1.165229 0.4504657 0.0000000
984 SARS-CoV-2 modulates host translation machinery 47 0.0000000 -0.7326514 -0.1199920 -0.4830142 -0.4257719 -0.6168887 0.0000000 0.1787794 0.0000000 0.0000002 0.0000000 1.160307 0.2318581 0.0000000
971 Ribosomal scanning and start codon recognition 57 0.0000000 -0.6334491 0.0824174 -0.5935089 -0.5536607 -0.5221832 0.0000000 0.2491757 0.0000000 0.0000000 0.0000000 1.157376 0.2972782 0.0000000
946 Regulation of ornithine decarboxylase (ODC) 49 0.0000000 -0.6932697 -0.6282254 0.1665645 0.2804466 -0.5844761 0.0000000 0.0000000 0.0267969 0.0002946 0.0000000 1.150346 0.4737044 0.0000000
1221 Translation initiation complex formation 57 0.0000000 -0.6168285 0.0867883 -0.5945860 -0.5555520 -0.5072511 0.0000000 0.2265856 0.0000000 0.0000000 0.0000000 1.143453 0.2960213 0.0000000
316 Erythropoietin activates Phosphoinositide-3-kinase (PI3K) 11 0.0025018 0.5834663 0.6385542 -0.3321371 -0.4501663 0.4952683 0.0005375 0.0002229 0.0527631 0.0080554 0.0031994 1.142996 0.5318835 0.0101625
621 Mitochondrial translation 94 0.0000000 -0.6589061 -0.6543841 0.1136986 0.2357808 -0.6049635 0.0000000 0.0000000 0.0324285 0.0000237 0.0000000 1.138805 0.4485159 0.0000000
196 Cooperation of Prefoldin and TriC/CCT in actin and tubulin folding 25 0.0000000 -0.6988509 -0.4677521 -0.3461136 -0.2128252 -0.6275953 0.0000000 0.0000336 0.0028535 0.0692309 0.0000000 1.125233 0.1991851 0.0000000
599 Metabolism of polyamines 54 0.0000000 -0.6481897 -0.6404774 0.2069916 0.3017422 -0.5471073 0.0000000 0.0000000 0.0042361 0.0000448 0.0000000 1.124091 0.4773296 0.0000000
44 Activation of the mRNA upon binding of the cap-binding complex and eIFs, and subsequent binding to 43S 58 0.0000000 -0.6214857 0.0688256 -0.5709758 -0.5314124 -0.5128567 0.0000000 0.3250797 0.0000000 0.0000000 0.0000000 1.123573 0.2839250 0.0000000
692 Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC) 113 0.0000000 -0.6314519 0.0988086 -0.5483800 -0.5102736 -0.5363730 0.0000000 0.0550157 0.0000000 0.0000000 0.0000000 1.121290 0.2966035 0.0000000
694 Nonsense-Mediated Decay (NMD) 113 0.0000000 -0.6314519 0.0988086 -0.5483800 -0.5102736 -0.5363730 0.0000000 0.0550157 0.0000000 0.0000000 0.0000000 1.121290 0.2966035 0.0000000
411 Generation of second messenger molecules 18 0.0000000 -0.3519137 0.3853261 -0.6286879 -0.6146235 -0.4491742 0.0117011 0.0041531 0.0000020 0.0000034 0.0011145 1.116731 0.4173076 0.0000000
934 Regulation of activated PAK-2p34 by proteasome mediated degradation 48 0.0000000 -0.6959441 -0.5737251 0.0784878 0.1863686 -0.5971654 0.0000000 0.0000000 0.2751249 0.0159677 0.0000000 1.100454 0.4176518 0.0000000
1115 Somitogenesis 46 0.0000000 -0.6688323 -0.5853209 0.1221077 0.2339324 -0.5900579 0.0000000 0.0000000 0.1091570 0.0032842 0.0000000 1.098972 0.4372670 0.0000000
787 Prefoldin mediated transfer of substrate to CCT/TriC 24 0.0000000 -0.6875960 -0.4495250 -0.3253785 -0.1877325 -0.6133934 0.0000000 0.0000910 0.0061292 0.1183004 0.0000001 1.091892 0.2045893 0.0000000
328 FBXL7 down-regulates AURKA during mitotic entry and in early mitosis 52 0.0000000 -0.6816868 -0.5846148 0.0856324 0.1961995 -0.5785338 0.0000000 0.0000000 0.2192121 0.0084083 0.0000000 1.089495 0.4178528 0.0000000
675 Negative regulation of NOTCH4 signaling 53 0.0000000 -0.6773215 -0.5766735 0.0547688 0.1754186 -0.5949953 0.0000000 0.0000000 0.4035726 0.0169357 0.0000000 1.085867 0.4046612 0.0000000
1264 Vif-mediated degradation of APOBEC3G 53 0.0000000 -0.6649836 -0.5912714 0.0828122 0.2034176 -0.5642545 0.0000000 0.0000000 0.2288520 0.0058710 0.0000000 1.076302 0.4146149 0.0000000
1254 Ubiquitin-dependent degradation of Cyclin D 50 0.0000000 -0.6733679 -0.5633933 0.1023407 0.2047238 -0.5768484 0.0000000 0.0000000 0.1562229 0.0070869 0.0000000 1.075163 0.4189392 0.0000000
1271 Vpu mediated degradation of CD4 51 0.0000000 -0.6824801 -0.5524564 0.0698858 0.1806959 -0.5855744 0.0000000 0.0000000 0.4378497 0.0287473 0.0000000 1.073042 0.4057351 0.0000000
1161 TP53 Regulates Transcription of Death Receptors and Ligands 10 0.0017944 0.2851102 -0.2311372 0.6656098 0.6384848 0.3908374 0.2185210 0.0775232 0.0002562 0.0004821 0.0660997 1.066848 0.3627780 0.0075958
199 Cross-presentation of soluble exogenous antigens (endosomes) 46 0.0000000 -0.6528929 -0.5698806 0.1262175 0.2205876 -0.5618807 0.0000000 0.0000000 0.0984802 0.0054885 0.0000000 1.063641 0.4236303 0.0000000
401 Gap junction assembly 10 0.0000451 -0.5248444 -0.1120360 -0.6347736 -0.5301865 -0.3979925 0.0028763 0.3188506 0.0008072 0.0046578 0.0160653 1.063225 0.2016147 0.0002655
1240 Transport of connexons to the plasma membrane 10 0.0000451 -0.5248444 -0.1120360 -0.6347736 -0.5301865 -0.3979925 0.0028763 0.3188506 0.0008072 0.0046578 0.0160653 1.063225 0.2016147 0.0002655
1253 Ubiquitin Mediated Degradation of Phosphorylated Cdc25A 49 0.0000000 -0.6687139 -0.5619567 0.0649142 0.1810283 -0.5745787 0.0000000 0.0000000 0.3517092 0.0179208 0.0000000 1.063061 0.4011894 0.0000000
1293 p53-Independent DNA Damage Response 49 0.0000000 -0.6687139 -0.5619567 0.0649142 0.1810283 -0.5745787 0.0000000 0.0000000 0.3517092 0.0179208 0.0000000 1.063061 0.4011894 0.0000000
1294 p53-Independent G1/S DNA damage checkpoint 49 0.0000000 -0.6687139 -0.5619567 0.0649142 0.1810283 -0.5745787 0.0000000 0.0000000 0.3517092 0.0179208 0.0000000 1.063061 0.4011894 0.0000000
961 Respiratory electron transport, ATP synthesis by chemiosmotic coupling, and heat production by uncoupling proteins. 123 0.0000000 -0.6324767 -0.5142019 -0.2122707 -0.1221275 -0.6350041 0.0000000 0.0000000 0.0000153 0.0123202 0.0000000 1.061901 0.2408708 0.0000000 
mypalette <- colorRampPalette(c("blue","white","red"))(n=25)
mmx <- head(mresm$enrichment_result,50)[,c(1,4:8)]
rownames(mmx) <- mmx[,1]
mmx[,1]=NULL
heatmap.2(as.matrix(mmx),scale="none",trace="none",
  margin=c(8,22), col=mypalette,cexCol=0.8)

# focus on apoptosis
mresm$enrichment_result[grep("apop",mresm$enrichment_result$set),] %>%
  kbl(caption="apoptosis pathways in productive compared to mock-infected") %>%
  kable_paper("hover",full_width=FALSE)
apoptosis pathways in productive compared to mock-infected
set setSize pMANOVA s.lat_v_by s.lat_v_prod s.prod_v_all s.prod_v_mock s.lat_v_mock p.lat_v_by p.lat_v_prod p.prod_v_all p.prod_v_mock p.lat_v_mock s.dist SD p.adjustMANOVA
136 Caspase activation via extrinsic apoptotic signalling pathway 21 0.1936749 0.1893086 0.1788480 0.0717101 0.0632074 0.2823470 0.1881514 0.2448373 0.5154099 0.5721466 0.0394881 0.3958304 0.0912475 0.3194536
359 Formation of apoptosome 10 0.3912958 0.1036601 -0.0072090 0.1923665 0.2183987 0.0330072 0.5202848 0.9954584 0.2504322 0.1991188 0.8061831 0.3107889 0.0977515 0.5097758
949 Regulation of the apoptosome activity 10 0.3912958 0.1036601 -0.0072090 0.1923665 0.2183987 0.0330072 0.5202848 0.9954584 0.2504322 0.1991188 0.8061831 0.3107889 0.0977515 0.5097758
205 Cytochrome c-mediated apoptotic response 12 0.1851160 0.1642874 0.0477389 0.1791861 0.1590466 0.0475600 0.2837528 0.7472528 0.2394247 0.3002126 0.7209370 0.2982193 0.0660630 0.3084481
1165 TP53 regulates transcription of several additional cell death genes whose specific roles in p53-dependent apoptosis remain uncertain 11 0.9792885 0.1221296 0.0923965 0.0887583 0.0719127 0.1491150 0.7343021 0.8847747 0.5772069 0.6532901 0.6233959 0.2423582 0.0306467 0.9822924
1153 TNFR1-induced proapoptotic signaling 23 0.0706112 0.0620582 0.0514603 0.0993134 0.0407975 0.1738252 0.5316731 0.6335161 0.3490543 0.6802378 0.1196193 0.2196410 0.0540802 0.1547409 
if (!file.exists("multi_mitch.html")) {
  mitch_report(res=mres5,outfile="multi_mitch.html",overwrite=TRUE)
}

# custom apop list
am <- mm[which(rownames(mm) %in% apop),,drop=FALSE]
am
##             lat_v_by   lat_v_prod   prod_v_all  prod_v_mock  lat_v_mock
## BAX     -30.76251252 -11.92870938  -9.03875061  -5.32460664 -21.3976648
## PTEN     21.32051053   4.92530925  14.92296985  11.27138605  16.3910356
## CYCS    -19.97075249 -31.55237715   9.36009013   4.74623047 -14.0007153
## ST6GAL1  19.92884203  13.53052980  -0.68162839  -1.37049203  16.4314650
## BCL2L11   8.61209030  10.56552806  -3.29009416  -3.82304227   3.5419553
## MEF2C     7.76316507   0.19693914  14.10634431   6.56896458   3.7846933
## BCL2      6.86792163   0.02468332  14.50872197   8.12345445   4.0604840
## CASP9     6.83892720   4.15765685  -0.26071214  -0.33859942   5.0047795
## CASP3     6.44055866   4.08280387  -0.41067524  -4.04225029   0.5560729
## CASP7     5.63773348   0.68857954   6.22697375   4.66258312   4.4844346
## CASP8     5.56686954   1.13464266   3.77896500  -1.85377861   3.0321282
## HRK       5.02172480   3.18041204  -0.12479074  -0.09151097   3.3197164
## SIRT1     4.48002582   1.20325251  -1.78976622  -0.64886280   2.1062324
## BCL2L2    4.14001595   0.25370945   5.28215310   2.47850099   2.0388448
## MCL1      3.77807893   4.27711111  -1.47518417  -2.22019479   1.4901105
## SELENOS  -3.71681388  -3.45887483  -0.03685491  -0.46553562  -3.8725593
## BCL2L1    3.49356694   0.25803111   5.49833946   3.58884604   2.4349476
## BID      -3.46668487  -9.77445785   9.38496521   8.71811788  -1.3141669
## BAD      -2.92209762 -10.94394139  11.99443548   8.90529930  -1.5041995
## BCL2L13   1.55612052  -1.53199007  12.43608061   7.85896810   1.1259759
## BAK1     -1.17457054  -1.67359684  -1.09534842   1.59402864  -0.3351071
## CCR5     -1.13990223   0.85072500  -0.03688655  -0.17021668  -0.5885146
## PLEKHO2  -0.87675679   1.33181829 -11.74504340 -16.28613144  -3.9381234
## BBC3      0.62738833  -0.39134438   3.96761906   4.09279560   1.2129514
## CDKN2A   -0.02386415   1.11444415  -2.58024064  -1.32598799  -0.1290172
heatmap.2(as.matrix(am),scale="none",trace="none",
  margin=c(8,22), col=mypalette,cexCol=0.8)

aresm <- mitch_calc(mm,genesets=list("custom_apop"=apop))
## Note: When prioritising by significance (ie: small 
##             p-values), large effect sizes might be missed.
t(aresm$enrichment_result)
##                1            
## set            "custom_apop"
## setSize        "25"         
## pMANOVA        "0.7283974"  
## s.lat_v_by     "-0.02114937"
## s.lat_v_prod   "-0.04584704"
## s.prod_v_all   "0.04961027" 
## s.prod_v_mock  "-0.01147897"
## s.lat_v_mock   "-0.08558379"
## p.lat_v_by     "0.9394897"  
## p.lat_v_prod   "0.7322448"  
## p.prod_v_all   "0.5977598"  
## p.prod_v_mock  "0.9672335"  
## p.lat_v_mock   "0.5144553"  
## s.dist         "0.1116547"  
## SD             "0.0496162"  
## p.adjustMANOVA "0.7283974"

Session information

For reproducibility

sessionInfo()
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] pkgload_1.3.2.1  GGally_2.1.2     ggplot2_3.4.3    reshape2_1.4.4  
##  [5] beeswarm_0.4.0   gtools_3.9.4     tibble_3.2.1     dplyr_1.1.3     
##  [9] echarts4r_0.4.5  kableExtra_1.3.4 gplots_3.1.3     mitch_1.12.0    
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.4       xfun_0.40          bslib_0.5.1        htmlwidgets_1.6.2 
##  [5] caTools_1.18.2     vctrs_0.6.3        tools_4.3.1        bitops_1.0-7      
##  [9] generics_0.1.3     parallel_4.3.1     fansi_1.0.4        highr_0.10        
## [13] pkgconfig_2.0.3    KernSmooth_2.23-22 RColorBrewer_1.1-3 webshot_0.5.5     
## [17] lifecycle_1.0.3    compiler_4.3.1     stringr_1.5.0      munsell_0.5.0     
## [21] httpuv_1.6.11      htmltools_0.5.6    sass_0.4.7         yaml_2.3.7        
## [25] later_1.3.1        pillar_1.9.0       jquerylib_0.1.4    MASS_7.3-60       
## [29] ellipsis_0.3.2     cachem_1.0.8       mime_0.12          tidyselect_1.2.0  
## [33] rvest_1.0.3        digest_0.6.33      stringi_1.7.12     fastmap_1.1.1     
## [37] grid_4.3.1         colorspace_2.1-0   cli_3.6.1          magrittr_2.0.3    
## [41] utf8_1.2.3         withr_2.5.0        scales_1.2.1       promises_1.2.1    
## [45] rmarkdown_2.25     httr_1.4.7         gridExtra_2.3      shiny_1.7.5       
## [49] evaluate_0.21      knitr_1.44         viridisLite_0.4.2  rlang_1.1.1       
## [53] Rcpp_1.0.11        xtable_1.8-4       glue_1.6.2         xml2_1.3.5        
## [57] svglite_2.1.1      rstudioapi_0.15.0  reshape_0.8.9      jsonlite_1.8.7    
## [61] R6_2.5.1           plyr_1.8.8         systemfonts_1.0.4