Lung-brain axis

Bioinformatics methods

Fastqc (v0.11.9) was used to inspect sequence quality[1].

The mouse transcriptome was downloaded from GENCODE version 28[2].

Skewer (v0.2.2) was used to trim low quality bases (qual<20) from the 3’ end of the read[3].

Kallisto (0.46.1) was used to map RNA-seq reads to the transcriptome [4].

Multiqc was used to tabulate sequence quality, trimming and mapping statistics [5].

Data were read into R v4.1.2 and duplicate lane data were aggregated, and transcript level counts were aggregated to gene level counts.

Genes with an average of less than 10 reads across samples were excluded from downstream analysis.

DESeq (1.32.0) was used with default settings to assay differential expression between control and treatment groups for all tissues [6].

Pathway analysis was performed with reactome gene sets obtained from MSigDB and converted to mouse gene identifiers with the msigdbr package (performed on 16-02-2022) [7,8,9].

Differential pathway analysis was performed with the “mitch” bioconductor package [10].

Genes and pathways with false discovery rate (FDR)<0.05 were considered significant.

suppressPackageStartupMessages({
    library("zoo")
    library("tidyverse")
    library("reshape2")
    library("DESeq2")
    library("gplots")
    library("fgsea")
    library("MASS")
    library("mitch")
    library("eulerr")
    library("limma")
    library("topconfects")
    library("kableExtra")
    library("vioplot")
    library("beeswarm")
})

Import read counts

Importing RNA-seq data

tmp <- read.table("3col.tsv.gz",header=F)
x <- as.matrix(acast(tmp, V2~V1, value.var="V3", fun.aggregate = sum))
x <- as.data.frame(x)
accession <- sapply((strsplit(rownames(x),"\\|")),"[[",2)
symbol<-sapply((strsplit(rownames(x),"\\|")),"[[",6)
x$geneid <- paste(accession,symbol)
xx <- aggregate(. ~ geneid,x,sum)
rownames(xx) <- xx$geneid
xx$geneid = NULL
xx <- round(xx)
xx <- xx[,which(colnames(xx)!="test")]
xx[1:6,1:6]

##                             EA11-1_S4_L001 EA11-1_S4_L002 EA11-2_S28_L001
## ENSMUSG00000000001.5 Gnai3             288            331             302
## ENSMUSG00000000003.16 Pbsn               0              0               0
## ENSMUSG00000000028.16 Cdc45             23             12              15
## ENSMUSG00000000031.17 H19                8             13               2
## ENSMUSG00000000037.18 Scml2             11              9               5
## ENSMUSG00000000049.12 Apoh               2              3               2
##                             EA11-2_S28_L002 EA11-3_S52_L001 EA11-3_S52_L002
## ENSMUSG00000000001.5 Gnai3              315             270             276
## ENSMUSG00000000003.16 Pbsn                0               0               0
## ENSMUSG00000000028.16 Cdc45              23              32              15
## ENSMUSG00000000031.17 H19                 4               3               5
## ENSMUSG00000000037.18 Scml2               3               5               9
## ENSMUSG00000000049.12 Apoh                0               4               3

dim(xx)

## [1] 55357   208

Fix the sample names.

They are duplicated for lane 1 and 2, which I will aggregate.

txx <- as.data.frame(t(xx))
txx$label <- sapply(strsplit(rownames(txx),"_"),"[[",1)
txx[1:3,c(1:4,ncol(txx))]

##                 ENSMUSG00000000001.5 Gnai3 ENSMUSG00000000003.16 Pbsn
## EA11-1_S4_L001                         288                          0
## EA11-1_S4_L002                         331                          0
## EA11-2_S28_L001                        302                          0
##                 ENSMUSG00000000028.16 Cdc45 ENSMUSG00000000031.17 H19  label
## EA11-1_S4_L001                           23                         8 EA11-1
## EA11-1_S4_L002                           12                        13 EA11-1
## EA11-2_S28_L001                          15                         2 EA11-2

txx2 <- aggregate(. ~ label,txx,sum)
txx2[1:3,1:4]

##    label ENSMUSG00000000001.5 Gnai3 ENSMUSG00000000003.16 Pbsn
## 1 EA11-1                        619                          0
## 2 EA11-2                        617                          0
## 3 EA11-3                        546                          0
##   ENSMUSG00000000028.16 Cdc45
## 1                          35
## 2                          38
## 3                          47

rownames(txx2) <- txx2[,1]
txx2[,1] = NULL
xx <- as.data.frame(t(txx2))
xx[1:4,1:5]

##                             EA11-1 EA11-2 EA11-3 EA11-4 EA11-5
## ENSMUSG00000000001.5 Gnai3     619    617    546    577    622
## ENSMUSG00000000003.16 Pbsn       0      0      0      0      0
## ENSMUSG00000000028.16 Cdc45     35     38     47     30     54
## ENSMUSG00000000031.17 H19       21      6      8     32     51

write.table(xx,file="lungbrain_counts.tsv",sep="\t",quote=FALSE)
rxx <- xx/colSums(xx) *1e6
rxx[1:4,1:5]

##                                EA11-1     EA11-2     EA11-3    EA11-4    EA11-5
## ENSMUSG00000000001.5 Gnai3  32.468943 31.6977874 30.2891348 26.133490 29.561201
## ENSMUSG00000000003.16 Pbsn   0.000000  0.0000000  0.0000000  0.000000  0.000000
## ENSMUSG00000000028.16 Cdc45  1.932574  2.0954246  2.6917763  1.428374  2.851489
## ENSMUSG00000000031.17 H19    1.181065  0.2705589  0.3625887  1.480344  2.481030

Samplesheet.

Need to delete “EA13-4”.

ss <- read.table("lungbrain_samplesheet.tsv",header=TRUE)

rownames(ss) <- gsub('\\.','-',paste("EA",as.character(ss$OriginalID),sep=""))
ss <- ss[order(rownames(ss)),]

rownames(ss) == colnames(rxx)

##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

ss %>% kbl(caption = "Sample sheet") %>%
  kable_paper("hover", full_width = F)

Sample sheet
	Region	Treatment	SampleNum	OriginalID	InternalID	PercentReads	NumReads	Gbases
EA11-1	hyp	OVA	4	11.1	CSEA004	0.5566	40914716	4.091472
EA11-2	amy	OVA	28	11.2	CSEA028	0.6315	46420487	4.642049
EA11-3	hip	OVA	52	11.3	CSEA052	0.5436	39959108	3.995911
EA11-4	pag	OVA	76	11.4	CSEA076	0.5279	38805028	3.880503
EA11-5	NI	OVA	32	11.5	CSEA032	0.6059	44538674	4.453867
EA12-1	hyp	OVA	5	12.1	CSEA005	0.6140	45134092	4.513409
EA12-2	amy	OVA	29	12.2	CSEA029	0.6140	45134092	4.513409
EA12-3	hip	OVA	53	12.3	CSEA053	0.4697	34526845	3.452684
EA12-4	pag	OVA	77	12.4	CSEA077	0.6622	48677191	4.867719
EA12-5	NI	OVA	40	12.5	CSEA040	0.6281	46170558	4.617056
EA13-1	hyp	Control	20	13.1	CSEA020	0.6535	48037669	4.803767
EA13-2	amy	Control	43	13.2	CSEA043	0.6126	45031180	4.503118
EA13-3	hip	Control	68	13.3	CSEA068	0.6158	45266406	4.526641
EA13-4	pag	Control	92	13.4	CSEA092	0.5902	43384594	4.338459
EA13-5	NI	Control	100	13.5	CSEA100	0.5576	40988224	4.098822
EA14-1	hyp	Control	21	14.1	CSEA021	0.6006	44149080	4.414908
EA14-2	amy	Control	44	14.2	CSEA044	0.6354	46707169	4.670717
EA14-4	pag	Control	93	14.4	CSEA093	0.6229	45788315	4.578832
EA14-5	NI	Control	101	14.5	CSEA101	0.4591	33747657	3.374766
EA15-3	hip	Control	69	15.3	CSEA069	0.6063	44568078	4.456808
EA16-2	amy	Control	45	16.2	CSEA045	0.5070	37268704	3.726870
EA19-1	hyp	Control	22	19.1	CSEA022	0.5334	39209323	3.920932
EA19-2	amy	Control	46	19.2	CSEA046	0.5431	39922354	3.992235
EA19-3	hip	Control	70	19.3	CSEA070	0.4972	36548323	3.654832
EA19-4	pag	Control	94	19.4	CSEA094	0.6218	45707456	4.570746
EA19-5	NI	Control	102	19.5	CSEA102	0.5271	38746221	3.874622
EA20-1	hyp	Control	23	20.1	CSEA023	0.5896	43340489	4.334049
EA20-2	amy	Control	47	20.2	CSEA047	0.5473	40231088	4.023109
EA20-3	hip	Control	71	20.3	CSEA071	0.5420	39841494	3.984149
EA20-4	pag	Control	95	20.4	CSEA095	0.6251	45950034	4.595003
EA20-5	NI	Control	103	20.5	CSEA103	0.5385	39584215	3.958422
EA22-5	NI	LPS	56	22.5	CSEA056	0.5462	40150229	4.015023
EA23-5	NI	LPS	64	23.5	CSEA064	0.7148	52543727	5.254373
EA24-1	hyp	LPS	9	24.1	CSEA009	0.6135	45097337	4.509734
EA24-2	amy	LPS	33	24.2	CSEA033	0.6361	46758625	4.675862
EA24-3	hip	LPS	57	24.3	CSEA057	0.5850	43002351	4.300235
EA24-4	pag	LPS	81	24.4	CSEA081	0.5860	43075859	4.307586
EA24-5	NI	LPS	72	24.5	CSEA072	0.5553	40819155	4.081915
EA25-1	hyp	LPS	10	25.1	CSEA010	0.6531	48008266	4.800827
EA25-2	amy	LPS	34	25.2	CSEA034	0.5554	40826506	4.082651
EA25-3	hip	LPS	58	25.3	CSEA058	0.5951	43744785	4.374478
EA25-4	pag	LPS	82	25.4	CSEA082	0.5224	38400732	3.840073
EA25-5	NI	LPS	80	25.5	CSEA080	0.5595	41127890	4.112789
EA26-1	hyp	LPS	11	26.1	CSEA011	0.5547	40775050	4.077505
EA26-2	amy	LPS	35	26.2	CSEA035	0.5051	37129038	3.712904
EA26-3	hip	LPS	59	26.3	CSEA059	0.5320	39106412	3.910641
EA26-4	pag	LPS	83	26.4	CSEA083	0.5903	43391945	4.339195
EA26-5	NI	LPS	88	26.5	CSEA088	0.5439	39981160	3.998116
EA28-1	hyp	LPS	12	28.1	CSEA012	0.5116	37606842	3.760684
EA28-2	amy	LPS	36	28.2	CSEA036	0.6100	44840058	4.484006
EA28-3	hip	LPS	60	28.3	CSEA060	0.5066	37239301	3.723930
EA28-4	pag	LPS	84	28.4	CSEA084	0.6240	45869174	4.586917
EA29-1	hyp	LPS	13	29.1	CSEA013	0.5781	42495144	4.249514
EA29-2	amy	LPS	37	29.2	CSEA037	0.5550	40797102	4.079710
EA29-3	hip	LPS	61	29.3	CSEA061	0.6357	46729221	4.672922
EA29-4	pag	LPS	85	29.4	CSEA085	0.5473	40231088	4.023109
EA29-5	NI	LPS	96	29.5	CSEA096	0.6012	44193185	4.419319
EA3-1	hyp	OVA	1	3.1	CSEA001	0.7106	52234993	5.223499
EA3-2	amy	OVA	25	3.2	CSEA025	0.5306	39003500	3.900350
EA3-3	hip	OVA	49	3.3	CSEA049	0.6483	47655426	4.765543
EA3-4	pag	OVA	73	3.4	CSEA073	0.5320	39106412	3.910641
EA3-5	NI	OVA	8	3.5	CSEA008	0.6426	47236429	4.723643
EA30-1	hyp	LPS	14	30.1	CSEA014	0.5787	42539249	4.253925
EA30-2	amy	LPS	38	30.2	CSEA038	0.5688	41811517	4.181152
EA30-3	hip	LPS	62	30.3	CSEA062	0.6359	46743923	4.674392
EA30-4	pag	LPS	86	30.4	CSEA086	0.6063	44568078	4.456808
EA30-5	NI	LPS	104	30.5	CSEA104	0.5692	41840920	4.184092
EA31-1	hyp	LPS	15	31.1	CSEA015	0.5738	42179058	4.217906
EA31-2	amy	LPS	39	31.2	CSEA039	0.5634	41414572	4.141457
EA31-3	hip	LPS	63	31.3	CSEA063	0.6144	45163495	4.516349
EA31-4	pag	LPS	87	31.4	CSEA087	0.5989	44024116	4.402412
EA34-1	hyp	OVA	6	34.1	CSEA006	0.6620	48662490	4.866249
EA34-2	amy	OVA	30	34.2	CSEA030	0.5227	38422784	3.842278
EA34-3	hip	OVA	54	34.3	CSEA054	0.6550	48147932	4.814793
EA34-4	pag	OVA	78	34.4	CSEA078	0.6248	45927981	4.592798
EA35-1	hyp	OVA	7	35.1	CSEA007	0.5053	37143740	3.714374
EA35-2	amy	OVA	31	35.2	CSEA031	0.5840	42928843	4.292884
EA35-3	hip	OVA	55	35.3	CSEA055	0.5617	41289608	4.128961
EA35-4	pag	OVA	79	35.4	CSEA079	0.6147	45185547	4.518555
EA35-5	NI	OVA	48	35.5	CSEA048	0.6017	44229940	4.422994
EA4-1	hyp	OVA	2	4.1	CSEA002	0.5883	43244928	4.324493
EA4-2	amy	OVA	26	4.2	CSEA026	0.5446	40032616	4.003262
EA4-3	hip	OVA	50	4.3	CSEA050	0.5873	43171420	4.317142
EA4-4	pag	OVA	74	4.4	CSEA074	0.5389	39613619	3.961362
EA4-5	NI	OVA	16	4.5	CSEA016	0.6581	48375807	4.837581
EA5-1	hyp	Control	17	5.1	CSEA017	0.5628	41370467	4.137047
EA5-2	amy	Control	41	5.2	CSEA041	0.5405	39731232	3.973123
EA5-3	hip	Control	65	5.3	CSEA065	0.6693	49199100	4.919910
EA5-4	pag	Control	89	5.4	CSEA089	0.5372	39488655	3.948865
EA5-5	NI	Control	97	5.5	CSEA097	0.6432	47280534	4.728053
EA6-1	hyp	Control	18	6.1	CSEA018	0.6405	47082061	4.708206
EA6-2	amy	Control	42	6.2	CSEA042	0.5692	41840920	4.184092
EA6-3	hip	Control	66	6.3	CSEA066	0.6563	48243492	4.824349
EA6-4	pag	Control	90	6.4	CSEA090	0.5525	40613332	4.061333
EA6-5	NI	Control	98	6.5	CSEA098	0.5549	40789751	4.078975
EA7-1	hyp	Control	19	7.1	CSEA019	0.5891	43303735	4.330374
EA7-3	hip	Control	67	7.3	CSEA067	0.6897	50698669	5.069867
EA7-4	pag	Control	91	7.4	CSEA091	0.5014	36857058	3.685706
EA7-5	NI	Control	99	7.5	CSEA099	0.6789	49904780	4.990478
EA9-1	hyp	OVA	3	9.1	CSEA003	0.5619	41304309	4.130431
EA9-2	amy	OVA	27	9.2	CSEA027	0.4965	36496867	3.649687
EA9-3	hip	OVA	51	9.3	CSEA051	0.4721	34703265	3.470326
EA9-4	pag	OVA	75	9.4	CSEA075	0.6358	46736572	4.673657
EA9-5	NI	OVA	24	9.5	CSEA024	0.5248	38577152	3.857715

QC analysis

Here I’ll look at a few different quality control measures.

The simplest is the number and proportion of reads that are assigned to genes.

par(mar=c(5,8,3,1))
barplot(colSums(xx),horiz=TRUE,las=1,xlab="num reads",col=ss$cols)

sums <- colSums(xx)
sums <- sums[order(sums)]
barplot(sums,horiz=TRUE,las=1,xlab="num reads",cex.names=0.8)
abline(v=15000000,col="red")

which(sums<15000000)

## EA9-5 
##     1

barplot(head(sums,6),horiz=TRUE,las=1,xlab="num reads",cex.names=1)
abline(v=15000000,col="red")

which(sums<15000000)

## EA9-5 
##     1

Now to look at the rRNA content.

rrna <- read.table("https://raw.githubusercontent.com/markziemann/lung_brain/main/rrna_res.tsv")
myrrna <- rrna$V2/10000
names(myrrna) <- rrna$V1
myrrna <- myrrna[order(myrrna)]

barplot(myrrna,horiz=TRUE,las=1,xlab="percent rRNA reads",cex.names=0.8)
abline(v=10,col="red")

barplot(tail(myrrna),horiz=TRUE,las=1,xlab="percent rRNA reads",cex.names=0.8)
abline(v=10,col="red")

myrrna[which(myrrna>10)]

##  EA29-5  EA23-5   EA3-1   EA9-5 
## 10.1406 15.5685 18.0218 25.9348

MDS plot for all samples

Multidimensional scaling plot to show the variation between all samples, very similar to PCA.

par(mar=c(5.1,4.1,4.1,2.1))

mds <- cmdscale(dist(t(xx)))

cols <- as.numeric(as.factor(ss$Region))+1
pchs <- as.numeric(factor(ss$Treatment))+14

plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n",pch=pchs, cex=4 ,col=cols )
text(mds, labels=rownames(mds) ,col="black")

legend("bottomleft", inset=.02, title="tissue",
   legend=unique(as.factor(ss$Region)) , fill=unique(cols),  cex=1.2)

legend("left", inset=.02, title="treatment",
   legend=unique(as.factor(ss$Treatment)) , pch=unique(pchs),  cex=1.2)

Correlation heatmap

heatmap.2(cor(xx),trace="n",main="Pearson correlation heatmap",margin=c(8,8),cexRow=0.7,cexCol=0.7)

heatmap.2(cor(xx),trace="n",main="Pearson correlation heatmap",margin=c(8,8),cexRow=0.7,cexCol=0.7)

Set up the different datasets for differential expression analysis

Here, I’ll give an example on how to separate the data matrix by tissue and then evaluate differential expression.

Don’t forget to remove poorly detected genes from the matrix with a threshold of 10 reads per sample on average.

There are 5 contrasts to set up, one for each tissue.

The separate sample sheets are called s1, s2, etc.

The separate counts tables are called x1, x2, etc.

I will begin with hypothalamus and leave the rest to Craig’s team.

# hyp 1
ss1 <- ss[which(ss$Region=="hyp"),]
xx1 <- xx[which(colnames(xx) %in% rownames(ss1))]
xx1 <- xx1[which(rowMeans(xx1)>=10),]
rpm1 <- xx1/colSums(xx1) *1e6
ss1_l <- ss[which(ss$Region=="hyp" & ss$Treatment != "OVA"),]
xx1_l <- xx[which(colnames(xx) %in% rownames(ss1_l))]
xx1_l <- xx1_l[which(rowMeans(xx1_l)>=10),]
rpm1_l <- xx1_l/colSums(xx1_l) *1e6
ss1_o <- ss[which(ss$Region=="hyp" & ss$Treatment != "LPS"),]
xx1_o <- xx[which(colnames(xx) %in% rownames(ss1_o))]
xx1_o <- xx1_o[which(rowMeans(xx1_o)>=10),]
rpm1_o <- xx1_o/colSums(xx1_o) *1e6

# amy 2
ss2 <- ss[which(ss$Region=="amy"),]
xx2 <- xx[which(colnames(xx) %in% rownames(ss2))]
xx2 <- xx2[which(rowMeans(xx2)>=10),]
rpm2 <- xx2/colSums(xx2) *1e6
ss2_l <- ss[which(ss$Region=="amy" & ss$Treatment != "OVA"),]
xx2_l <- xx[which(colnames(xx) %in% rownames(ss2_l))]
xx2_l <- xx2_l[which(rowMeans(xx2_l)>=10),]
rpm2_l <- xx2_l/colSums(xx2_l) *1e6
ss2_o <- ss[which(ss$Region=="amy" & ss$Treatment != "LPS"),]
xx2_o <- xx[which(colnames(xx) %in% rownames(ss2_o))]
xx2_o <- xx2_o[which(rowMeans(xx2_o)>=10),]
rpm1_o <- xx2_o/colSums(xx2_o) *1e6

# hip 3
ss3 <- ss[which(ss$Region=="hip"),]
xx3 <- xx[which(colnames(xx) %in% rownames(ss3))]
xx3 <- xx3[which(rowMeans(xx3)>=10),]
rpm3 <- xx3/colSums(xx3) *1e6
ss3_l <- ss[which(ss$Region=="hip" & ss$Treatment != "OVA"),]
xx3_l <- xx[which(colnames(xx) %in% rownames(ss3_l))]
xx3_l <- xx3_l[which(rowMeans(xx3_l)>=10),]
rpm3_l <- xx3_l/colSums(xx3_l) *1e6
ss3_o <- ss[which(ss$Region=="hip" & ss$Treatment != "LPS"),]
xx3_o <- xx[which(colnames(xx) %in% rownames(ss3_o))]
xx3_o <- xx3_o[which(rowMeans(xx3_o)>=10),]
rpm3_o <- xx3_o/colSums(xx3_o) *1e6

# pag 4
ss4 <- ss[which(ss$Region=="pag"),]
xx4 <- xx[which(colnames(xx) %in% rownames(ss4))]
xx4 <- xx4[which(rowMeans(xx4)>=10),]
rpm4 <- xx4/colSums(xx4) *1e6
ss4_l <- ss[which(ss$Region=="pag" & ss$Treatment != "OVA"),]
xx4_l <- xx[which(colnames(xx) %in% rownames(ss4_l))]
xx4_l <- xx4_l[which(rowMeans(xx4_l)>=10),]
rpm4_l <- xx4_l/colSums(xx4_l) *1e6
ss4_o <- ss[which(ss$Region=="pag" & ss$Treatment != "LPS"),]
xx4_o <- xx[which(colnames(xx) %in% rownames(ss4_o))]
xx4_o <- xx4_o[which(rowMeans(xx4_o)>=10),]
rpm4_o <- xx4_o/colSums(xx4_o) *1e6

# ni 5
ss5 <- ss[which(ss$Region=="NI"),]
xx5 <- xx[which(colnames(xx) %in% rownames(ss5))]
xx5 <- xx5[which(rowMeans(xx5)>=10),]
rpm5 <- xx5/colSums(xx5) *1e6
ss5_l <- ss[which(ss$Region=="NI" & ss$Treatment != "OVA"),]
xx5_l <- xx[which(colnames(xx) %in% rownames(ss5_l))]
xx5_l <- xx5_l[which(rowMeans(xx5_l)>=10),]
rpm5_l <- xx5_l/colSums(xx5_l) *1e6
ss5_o <- ss[which(ss$Region=="NI" & ss$Treatment != "LPS"),]
xx5_o <- xx[which(colnames(xx) %in% rownames(ss5_o))]
xx5_o <- xx5_o[which(rowMeans(xx5_o)>=10),]
rpm5_o <- xx5_o/colSums(xx5_o) *1e6

message("dimensions of objects: n genes, n samples")

## dimensions of objects: n genes, n samples

l <- list("xx"=xx,"xx1"=xx1,"xx1 l"= xx1_l,"xx1 o"=xx1_o,
  "xx2"=xx2,"xx2 l"=xx2_l,"xx2 o"=xx2_o,
  "xx3"=xx3,"xx3 l"=xx3_l,"xx3 o"=xx3_o,
  "xx4"=xx4,"xx4 l"=xx4_l,"xx4 o"=xx5_o,
  "xx5"=xx5,"xx5 l"=xx5_l,"xx5 o"=xx5_o)

lapply(l,dim)

## $xx
## [1] 55357   104
## 
## $xx1
## [1] 17091    21
## 
## $`xx1 l`
## [1] 17095    14
## 
## $`xx1 o`
## [1] 17114    14
## 
## $xx2
## [1] 16681    21
## 
## $`xx2 l`
## [1] 16818    14
## 
## $`xx2 o`
## [1] 16624    14
## 
## $xx3
## [1] 16897    21
## 
## $`xx3 l`
## [1] 16927    14
## 
## $`xx3 o`
## [1] 16953    14
## 
## $xx4
## [1] 17570    21
## 
## $`xx4 l`
## [1] 17612    14
## 
## $`xx4 o`
## [1] 17391    13
## 
## $xx5
## [1] 17474    20
## 
## $`xx5 l`
## [1] 17490    14
## 
## $`xx5 o`
## [1] 17391    13

Differential expression with DESeq2

For each brain region, there are two treatments, and so two contrasts are performed.

Before we do that, let’s look at the MDS plot of the hypothalamus data.

par(mar=c(5.1,4.1,4.1,2.1))

mds <- cmdscale(dist(t(xx1)))

cols <- as.numeric(factor(ss1$Treatment))+1

plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
  type = "p",bty="n", pch=19, cex=4 ,col=cols )
text(mds, labels=rownames(mds) ,col="black")

legend("bottomleft", inset=.02, title="treatment",
   legend=unique(as.factor(ss1$Treatment)) , col=unique(cols), pch=19,  cex=1.2)

Now run DESeq2 for all contrasts.

# 1 hyp
dds <- DESeqDataSetFromMatrix(countData = xx1_l , colData = ss1_l , design = ~ Treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 7 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge1_l <- dge
d1up_l <- rownames(subset(dge1_l,padj <= 0.05 & log2FoldChange > 0))
d1dn_l <- rownames(subset(dge1_l,padj <= 0.05 & log2FoldChange < 0))
write.table(dge1_l,file="dge1_l.tsv",quote=FALSE,sep="\t")

dds <- DESeqDataSetFromMatrix(countData = xx1_o , colData = ss1_o , design = ~ Treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 10 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge1_o <- dge
d1up_o <- rownames(subset(dge1_o,padj <= 0.05 & log2FoldChange > 0))
d1dn_o <- rownames(subset(dge1_o,padj <= 0.05 & log2FoldChange < 0))
write.table(dge1_o,file="dge1_o.tsv",quote=FALSE,sep="\t")

# 2 amy
dds <- DESeqDataSetFromMatrix(countData = xx2_l , colData = ss2_l , design = ~ Treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 5 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge2_l <- dge
d2up_l <- rownames(subset(dge2_l,padj <= 0.05 & log2FoldChange > 0))
d2dn_l <- rownames(subset(dge2_l,padj <= 0.05 & log2FoldChange < 0))
write.table(dge2_l,file="dge2_l.tsv",quote=FALSE,sep="\t")

dds <- DESeqDataSetFromMatrix(countData = xx2_o , colData = ss2_o , design = ~ Treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 19 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge2_o <- dge
d2up_o <- rownames(subset(dge2_o,padj <= 0.05 & log2FoldChange > 0))
d2dn_o <- rownames(subset(dge2_o,padj <= 0.05 & log2FoldChange < 0))
write.table(dge2_o,file="dge2_o.tsv",quote=FALSE,sep="\t")

# 3 hip
dds <- DESeqDataSetFromMatrix(countData = xx3_l , colData = ss3_l , design = ~ Treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 2 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge3_l <- dge
d3up_l <- rownames(subset(dge3_l,padj <= 0.05 & log2FoldChange > 0))
d3dn_l <- rownames(subset(dge3_l,padj <= 0.05 & log2FoldChange < 0))
write.table(dge3_l,file="dge3_l.tsv",quote=FALSE,sep="\t")

dds <- DESeqDataSetFromMatrix(countData = xx3_o , colData = ss3_o , design = ~ Treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 6 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge3_o <- dge
d3up_o <- rownames(subset(dge3_o,padj <= 0.05 & log2FoldChange > 0))
d3dn_o <- rownames(subset(dge3_o,padj <= 0.05 & log2FoldChange < 0))
write.table(dge3_o,file="dge3_o.tsv",quote=FALSE,sep="\t")

# 4 pag
dds <- DESeqDataSetFromMatrix(countData = xx4_l , colData = ss4_l , design = ~ Treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 6 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge4_l <- dge
d4up_l <- rownames(subset(dge4_l,padj <= 0.05 & log2FoldChange > 0))
d4dn_l <- rownames(subset(dge4_l,padj <= 0.05 & log2FoldChange < 0))
write.table(dge4_l,file="dge4_l.tsv",quote=FALSE,sep="\t")

dds <- DESeqDataSetFromMatrix(countData = xx4_o , colData = ss4_o , design = ~ Treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 4 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge4_o <- dge
d4up_o <- rownames(subset(dge4_o,padj <= 0.05 & log2FoldChange > 0))
d4dn_o <- rownames(subset(dge4_o,padj <= 0.05 & log2FoldChange < 0))
write.table(dge4_o,file="dge4_o.tsv",quote=FALSE,sep="\t")

# 5 NI
dds <- DESeqDataSetFromMatrix(countData = xx5_l , colData = ss5_l , design = ~ Treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 15 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge5_l <- dge
d5up_l <- rownames(subset(dge5_l,padj <= 0.05 & log2FoldChange > 0))
d5dn_l <- rownames(subset(dge5_l,padj <= 0.05 & log2FoldChange < 0))
write.table(dge5_l,file="dge5_l.tsv",quote=FALSE,sep="\t")

dds <- DESeqDataSetFromMatrix(countData = xx5_o , colData = ss5_o , design = ~ Treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

res <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 41 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
dge5_o <- dge
d5up_o <- rownames(subset(dge5_o,padj <= 0.05 & log2FoldChange > 0))
d5dn_o <- rownames(subset(dge5_o,padj <= 0.05 & log2FoldChange < 0))
write.table(dge5_o,file="dge5_o.tsv",quote=FALSE,sep="\t")

Here let’s look at some plots.

MA plot shows the average level and fold change of all detected genes. Volcano plot shows the fold change and the significance, as measured by -log(p-value). Significant genes are shown as red points.

There are heatmaps of the top ranked genes by p-value. Above the gene expression values there is a bar in orange/gray colours. Control is shown in orange and treatment in grey.

maplot <- function(de,contrast_name) {
  de <- de[which(!is.na(de$padj)),]
  sig <-subset(de, padj < 0.05 )
  up <-rownames(subset(de, padj < 0.05 & log2FoldChange > 0))
  dn <-rownames(subset(de, padj < 0.05 & log2FoldChange < 0))
  GENESUP <- length(up)
  GENESDN <- length(dn)
  DET=nrow(de)
  SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down", DET, "detected")
  ns <-subset(de, padj > 0.05 )
  plot(log2(de$baseMean),de$log2FoldChange,
       xlab="log2 basemean", ylab="log2 foldchange",
       pch=19, cex=0.5, col="dark gray",
       main=contrast_name, cex.main=1)
  points(log2(sig$baseMean),sig$log2FoldChange,
         pch=19, cex=0.5, col="red")
  mtext(SUBHEADER,cex = 1)
}

make_volcano <- function(de,name) {
    de <- de[which(!is.na(de$padj)),]
    de$pvalue[which(de$pvalue==0)] <- 1e-320
    sig <- subset(de,padj<0.05)
    N_SIG=nrow(sig)
    N_UP=nrow(subset(sig,log2FoldChange>0))
    N_DN=nrow(subset(sig,log2FoldChange<0))
    DET=nrow(de)
    HEADER=paste(N_SIG,"@5%FDR,", N_UP, "up", N_DN, "dn", DET, "detected")
    plot(de$log2FoldChange,-log10(de$pval),cex=0.5,pch=19,col="darkgray",
        main=name, xlab="log2 FC", ylab="-log10 pval")
    mtext(HEADER)
    grid()
    points(sig$log2FoldChange,-log10(sig$pval),cex=0.5,pch=19,col="red")
}

make_volcano2 <- function(de,name) {
    de <- de[which(!is.na(de$padj)),]
    de$pvalue[which(de$pvalue==0)] <- 1e-320
    sig <- subset(de,padj<0.05)
    N_SIG=nrow(sig)
    N_UP=nrow(subset(sig,log2FoldChange>0))
    N_DN=nrow(subset(sig,log2FoldChange<0))
    DET=nrow(de)
    HEADER=paste(N_SIG,"@5%FDR,", N_UP, "up", N_DN, "dn", DET, "detected")
    top <- head(sig,30)
    mylabels <- sapply(strsplit(rownames(top)," "),"[[",2)
    plot(de$log2FoldChange,-log10(de$pval),cex=0.5,pch=19,col="darkgray",
        main=name, xlab="log2 FC", ylab="-log10 pval")
    mtext(HEADER)
    grid()
    points(sig$log2FoldChange,-log10(sig$pval),cex=0.5,pch=19,col="red")
    text(top$log2FoldChange+0.2,-log10(top$pval),labels=mylabels, srt=35 ,cex=0.7)
}

make_heatmap <- function(de,name,myss,mx,n=30){
  colfunc <- colorRampPalette(c("blue", "white", "red"))
  csc <- myss$Treatment
  csc <- gsub("Control","orange",csc)
  csc <- gsub("LPS","gray",csc)
  csc <- gsub("OVA","yellow",csc)
  mxn <- mx/rowSums(mx)*1000000
  x <- mxn[which(rownames(mxn) %in% rownames(head(de,n))),]
  heatmap.2(as.matrix(x),trace="none",col=colfunc(25),scale="row", margins = c(10,15), cexRow=0.7,
    main=paste("Top ranked",n,"genes in",name) , ColSideColors = csc  )
  mtext("ctrl=orange, LPS=gray, OVA=yellow")
}

make_heatmap2 <- function(de,name,myss,mx,n=30){
  colfunc <- colorRampPalette(c("blue", "white", "red"))
  csc <- myss$Treatment
  csc <- gsub("Control","orange",csc)
  csc <- gsub("LPS","gray",csc)
  csc <- gsub("OVA","yellow",csc)
  mxn <- mx/rowSums(mx)*1000000
  x <- mxn[which(rownames(mxn) %in% rownames(head(de,n))),]
  rownames(x) <- sapply(strsplit(rownames(x)," "),"[[",2)
  heatmap.2(as.matrix(x),trace="none",col=colfunc(25),scale="row", margins = c(10,15), cexRow=0.7,
    main=paste("Top ranked",n,"genes in",name) ,  ColSideColors = csc   )
  mtext("ctrl=orange, LPS=gray, OVA=yellow")
}

mymds <- function(de,name,myss,mx) {
  mds <- cmdscale(dist(t(mx)))
  csc <-  myss$Treatment
  csc <- gsub("Control","orange",csc)
  csc <- gsub("LPS","gray",csc)
  csc <- gsub("OVA","yellow",csc)
  plot(mds, xlab="Coordinate 1", ylab="Coordinate 2", main = name ,
    type = "p",bty="n",pch=19, cex=4 ,col=csc )
  text(mds, labels=rownames(mds) ,col="black")
  legend("topright", inset=.02, title="treatment",
    legend=unique(as.factor(myss$Treatment)) , pch=19, col=unique(csc),  cex=1.4)
}

Results for hypothalamus

We show an MDS plot, MA plot, volcano plot, heatmap and table of top results for each contrast.

mymds(de=dge1_l,name="Cont1: Effect of LPS treatment in hypothalamus",myss=ss1_l,mx=xx1_l)

maplot(dge1_l,"Cont1: Effect of LPS treatment in hypothalamus")

make_volcano(dge1_l,"Cont1: Effect of LPS treatment in hypothalamus")

make_heatmap(de=dge1_l,name="Cont1: Effect of LPS treatment in hypothalamus",myss=ss1_l,mx=xx1_l,n=50)

make_heatmap2(de=dge1_l,name="Cont1: Effect of LPS treatment in hypothalamus",myss=ss1_l,mx=xx1_l,n=50)

dge1_l[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 1: Effect of LPS treatment in hypothalamus") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of LPS treatment in hypothalamus
	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000054252.19 Fgfr3	2892.62590	0.3328185	0.0609055	5.464504	0.00e+00	0.0006909
ENSMUSG00000022715.4 Tmem114	100.67950	0.8836136	0.1650244	5.354442	1.00e-07	0.0006909
ENSMUSG00000025498.16 Irf7	77.47802	-0.8571473	0.1692562	-5.064199	4.00e-07	0.0017092
ENSMUSG00000027875.13 Hmgcs2	180.91714	-0.6653837	0.1315622	-5.057562	4.00e-07	0.0017092
ENSMUSG00000028820.14 Sfpq	2569.85119	0.2618911	0.0527559	4.964201	7.00e-07	0.0022213
ENSMUSG00000041926.16 Rnpep	654.35913	0.2935820	0.0598784	4.902971	9.00e-07	0.0025330
ENSMUSG00000079363.8 Gbp4	72.55881	-1.1002191	0.2263856	-4.859934	1.20e-06	0.0027008
ENSMUSG00000018822.8 Sfrp5	659.14853	0.4883485	0.1033004	4.727461	2.30e-06	0.0041077
ENSMUSG00000026222.17 Sp100	55.19664	-0.7479415	0.1582800	-4.725433	2.30e-06	0.0041077
ENSMUSG00000028655.12 Mfsd2a	443.11911	0.5572282	0.1205845	4.621059	3.80e-06	0.0061467
ENSMUSG00000040253.16 Gbp7	103.19267	-0.6006230	0.1353212	-4.438500	9.10e-06	0.0132588
ENSMUSG00000004891.17 Nes	206.29657	0.4958473	0.1122069	4.419045	9.90e-06	0.0133010
ENSMUSG00000035578.17 Iqcg	381.27559	0.4608992	0.1047871	4.398437	1.09e-05	0.0135034
ENSMUSG00000108852.2 Gm44911	34.74705	0.9903775	0.2263409	4.375601	1.21e-05	0.0139263
ENSMUSG00000024907.8 Gal	1601.34229	0.3959929	0.0914229	4.331441	1.48e-05	0.0153554
ENSMUSG00000020717.20 Pecam1	340.34994	-0.3114213	0.0720066	-4.324902	1.53e-05	0.0153554
ENSMUSG00000051159.17 Cited1	338.53219	0.6483567	0.1506781	4.302925	1.69e-05	0.0159634
ENSMUSG00000051427.15 Ccdc157	914.73405	0.2973665	0.0699150	4.253256	2.11e-05	0.0188445
ENSMUSG00000028519.17 Dab1	925.88165	-0.2341632	0.0557089	-4.203335	2.63e-05	0.0222867
ENSMUSG00000064247.15 Plcxd1	838.55340	0.2157020	0.0515714	4.182586	2.88e-05	0.0232011

mymds(de=dge1_o,name="Cont1: Effect of OVA treatment in hypothalamus",myss=ss1_o,mx=xx1_o)

maplot(dge1_o,"Cont1: Effect of OVA treatment in hypothalamus")

make_volcano(dge1_o,"Cont1: Effect of OVA treatment in hypothalamus")

make_heatmap(de=dge1_o,name="Cont1: Effect of OVA treatment in hypothalamus",myss=ss1_o,mx=xx1_o,n=50)

make_heatmap2(de=dge1_o,name="Cont1: Effect of OVA treatment in hypothalamus",myss=ss1_o,mx=xx1_o,n=50)

dge1_o[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 1: Effect of treatment in hypothalamus") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 1: Effect of treatment in hypothalamus
	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000078880.11 Gm14308	92.882830	-0.9356691	0.2192963	-4.266689	0.0000198	0.3395346
ENSMUSG00000036915.18 Kirrel2	169.164941	0.8427760	0.2297223	3.668673	0.0002438	0.9998282
ENSMUSG00000044285.9 Ubb-ps	2336.340595	-0.3680306	0.1025092	-3.590221	0.0003304	0.9998282
ENSMUSG00000085328.3 Gm17131	141.467376	-0.6729148	0.1891178	-3.558178	0.0003734	0.9998282
ENSMUSG00000030413.8 Pglyrp1	106.131949	0.6503273	0.1833759	3.546416	0.0003905	0.9998282
ENSMUSG00000033502.16 Cdc14a	100.974101	0.4596276	0.1301010	3.532853	0.0004111	0.9998282
ENSMUSG00000108601.2 Gm44645	94.446899	0.6303080	0.1803630	3.494664	0.0004747	0.9998282
ENSMUSG00000075014.2 Gm10800	33.751721	3.2629734	0.9463915	3.447805	0.0005652	0.9998282
ENSMUSG00000054252.19 Fgfr3	2730.435705	0.2147599	0.0630824	3.404437	0.0006630	0.9998282
ENSMUSG00000059970.8 Hspa2	641.897572	0.3132088	0.0939489	3.333822	0.0008566	0.9998282
ENSMUSG00000092483.2 Gm20421	22.954083	0.7081938	0.2204315	3.212761	0.0013147	0.9998282
ENSMUSG00000021892.15 Sh3bp5	762.980552	-0.1810008	0.0565574	-3.200302	0.0013728	0.9998282
ENSMUSG00000041079.13 Rwdd2b	175.004412	0.2696634	0.0866861	3.110802	0.0018658	0.9998282
ENSMUSG00000044707.8 Ccnjl	31.491228	-0.6302744	0.2031683	-3.102228	0.0019207	0.9998282
ENSMUSG00000053279.9 Aldh1a1	339.445293	0.4826510	0.1581715	3.051441	0.0022775	0.9998282
ENSMUSG00000020123.7 Avpr1a	98.069219	0.4957189	0.1626906	3.047004	0.0023113	0.9998282
ENSMUSG00000020374.17 Rasgef1c	258.064242	-0.2698852	0.0889452	-3.034285	0.0024111	0.9998282
ENSMUSG00000091318.3 Gm5415	9.042075	1.4281893	0.4797236	2.977109	0.0029098	0.9998282
ENSMUSG00000021750.17 Fam107a	3483.985160	0.2995835	0.1007207	2.974398	0.0029356	0.9998282
ENSMUSG00000013089.16 Etv5	587.848238	0.2785646	0.0939914	2.963723	0.0030394	0.9998282

Results for amygdala

mymds(de=dge2_l,name="Cont2: Effect of LPS treatment in amygdala",myss=ss2_l,mx=xx2_l)

maplot(dge2_l,"Cont2: Effect of LPS treatment in amygdala")

make_volcano(dge2_l,"Cont2: Effect of LPS treatment in amygdala")

make_heatmap(de=dge2_l,name="Cont2: Effect of LPS treatment in amygdala",myss=ss2_l,mx=xx2_l,n=50)

make_heatmap2(de=dge2_l,name="Cont2: Effect of LPS treatment in amygdala",myss=ss2_l,mx=xx2_l,n=50)

dge2_l[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 2: Effect of LPS treatment in amygdala") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 2: Effect of LPS treatment in amygdala
	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000025498.16 Irf7	54.87672	-1.2294874	0.2232120	-5.508159	0.0000000	0.0006098
ENSMUSG00000079363.8 Gbp4	62.26653	-0.9118851	0.1718895	-5.305066	0.0000001	0.0006560
ENSMUSG00000036915.18 Kirrel2	254.72910	0.6575428	0.1241095	5.298088	0.0000001	0.0006560
ENSMUSG00000032902.2 Slc16a1	429.58840	-0.4275636	0.0835198	-5.119305	0.0000003	0.0011124
ENSMUSG00000018822.8 Sfrp5	179.17105	0.6336687	0.1241260	5.105042	0.0000003	0.0011124
ENSMUSG00000078921.4 Tgtp2	10.75353	-2.2904610	0.4554595	-5.028902	0.0000005	0.0013827
ENSMUSG00000027875.13 Hmgcs2	130.58206	-0.7567441	0.1547343	-4.890602	0.0000010	0.0024153
ENSMUSG00000007097.15 Atp1a2	13121.54535	0.2852928	0.0589693	4.837987	0.0000013	0.0027573
ENSMUSG00000027333.19 Smox	993.08881	0.3229665	0.0686244	4.706291	0.0000025	0.0047140
ENSMUSG00000004891.17 Nes	147.63914	0.5842943	0.1253589	4.660972	0.0000031	0.0052929
ENSMUSG00000102275.2 Gm37144	83.14632	-1.1552394	0.2515488	-4.592506	0.0000044	0.0066959
ENSMUSG00000089824.11 Rbm12	379.55172	-0.3020121	0.0696401	-4.336755	0.0000145	0.0202659
ENSMUSG00000060862.11 Zbtb40	209.30249	0.4768500	0.1107158	4.306972	0.0000166	0.0214112
ENSMUSG00000044712.16 Slc38a6	1014.61207	-0.4381471	0.1034282	-4.236245	0.0000227	0.0273039
ENSMUSG00000020680.16 Taf15	956.83934	0.3203512	0.0800847	4.000154	0.0000633	0.0709734
ENSMUSG00000034173.14 Zbed5	124.05394	-0.4521479	0.1157090	-3.907630	0.0000932	0.0979709
ENSMUSG00000006390.16 Elovl1	217.76918	0.3244737	0.0834422	3.888604	0.0001008	0.0997431
ENSMUSG00000052911.10 Lamb2	645.20070	0.3401842	0.0883340	3.851114	0.0001176	0.1098606
ENSMUSG00000031465.7 Angpt2	32.01754	-0.7723745	0.2036123	-3.793359	0.0001486	0.1315547
ENSMUSG00000100750.2 Gm29084	513.48571	-0.3339886	0.0899518	-3.712974	0.0002048	0.1722483

mymds(de=dge2_o,name="Cont2: Effect of OVA treatment in amygdala",myss=ss2_o,mx=xx2_o)

maplot(dge2_o,"Cont2: Effect of OVA treatment in amygdala")

make_volcano(dge2_o,"Cont2: Effect of OVA treatment in amygdala")

make_heatmap(de=dge2_o,name="Cont2: Effect of OVA treatment in amygdala",myss=ss2_o,mx=xx2_o,n=50)

make_heatmap2(de=dge2_o,name="Cont2: Effect of OVA treatment in amygdala",myss=ss2_o,mx=xx2_o,n=50)

dge2_o[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 2: Effect of treatment in amygdala") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 2: Effect of treatment in amygdala
	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000107676.2 Gm44178	185.4223	-0.8648437	0.1423153	-6.076956	0.0e+00	0.0000164
ENSMUSG00000025791.19 Pgm1	1679.6378	0.2939654	0.0506255	5.806664	0.0e+00	0.0000315
ENSMUSG00000030889.15 Vwa3a	101.9254	-0.7011680	0.1211103	-5.789499	0.0e+00	0.0000315
ENSMUSG00000059149.18 Mfsd4a	841.7077	-0.3390365	0.0621451	-5.455561	0.0e+00	0.0001636
ENSMUSG00000036915.18 Kirrel2	266.3206	0.8297410	0.1549163	5.356062	1.0e-07	0.0002019
ENSMUSG00000024782.8 Ak3	1076.8616	0.3287706	0.0615097	5.345017	1.0e-07	0.0002019
ENSMUSG00000018822.8 Sfrp5	171.4710	0.6007642	0.1168718	5.140370	3.0e-07	0.0005249
ENSMUSG00000032507.16 Fbxl2	973.7927	-0.3006599	0.0592058	-5.078214	4.0e-07	0.0006064
ENSMUSG00000053279.9 Aldh1a1	232.1522	0.5728101	0.1130798	5.065538	4.0e-07	0.0006064
ENSMUSG00000034413.15 Neurl1b	998.7317	-0.4392775	0.0892047	-4.924376	8.0e-07	0.0010672
ENSMUSG00000038633.6 Degs1	1216.6775	0.3005384	0.0611145	4.917627	9.0e-07	0.0010672
ENSMUSG00000030256.12 Bhlhe41	1281.0403	0.3288133	0.0684756	4.801908	1.6e-06	0.0017551
ENSMUSG00000030499.10 Kctd15	178.8979	-0.4463070	0.0942096	-4.737384	2.2e-06	0.0022317
ENSMUSG00000030209.15 Grin2b	4116.6194	-0.4866258	0.1043231	-4.664600	3.1e-06	0.0029599
ENSMUSG00000038738.16 Shank1	6687.6830	-0.3929989	0.0868324	-4.525946	6.0e-06	0.0053717
ENSMUSG00000006390.16 Elovl1	227.2895	0.5070433	0.1126133	4.502516	6.7e-06	0.0054351
ENSMUSG00000038260.11 Trpm4	233.7264	-0.4194385	0.0935316	-4.484457	7.3e-06	0.0054351
ENSMUSG00000034818.18 Celf5	4451.3635	-0.3143020	0.0701827	-4.478343	7.5e-06	0.0054351
ENSMUSG00000079018.11 Ly6c1	374.3433	0.4597376	0.1027940	4.472418	7.7e-06	0.0054351
ENSMUSG00000031513.9 Leprotl1	1221.1671	0.2195997	0.0492131	4.462220	8.1e-06	0.0054351

Results for hippocampus

mymds(de=dge3_l,name="Cont3: Effect of LPS treatment in hippocampus",myss=ss3_l,mx=xx3_l)

maplot(dge3_l,"Cont3: Effect of LPS treatment in hippocampus")

make_volcano(dge3_l,"Cont3: Effect of LPS treatment in hippocampus")

make_heatmap(de=dge3_l,name="Cont3: Effect of LPS treatment in hippocampus",myss=ss3_l,mx=xx3_l,n=50)

make_heatmap2(de=dge3_l,name="Cont3: Effect of LPS treatment in hippocampus",myss=ss3_l,mx=xx3_l,n=50)

dge3_l[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 3: Effect of LPS treatment in hippocampus") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 3: Effect of LPS treatment in hippocampus
	baseMean	log2FoldChange	lfcSE	stat
ENSMUSG00000100862.2 Gm10925	4975.55176	7.1868883	0.6335748	11.343393
ENSMUSG00000070713.6 Gm10282	208.77183	-1.0025619	0.0983370	-10.195161
ENSMUSG00000040447.16 Spns2	1554.46939	0.5308742	0.0527511	10.063750
ENSMUSG00000081603.2 Gm14681	107.33155	-1.1109529	0.1168663	-9.506187
ENSMUSG00000068882.14 Ssb	1815.27129	-0.5934983	0.0624491	-9.503708
ENSMUSG00000096006.2 Gm21596	253.23289	-0.7510200	0.0804176	-9.338997
ENSMUSG00000024261.7 Syt4	2094.97545	-0.5610059	0.0603557	-9.294993
ENSMUSG00000025104.14 Hdgfl3	2080.41389	-0.4435812	0.0489446	-9.062924
ENSMUSG00000017831.9 Rab5a	2217.40639	-0.4345348	0.0480838	-9.037037
ENSMUSG00000047879.18 Usp14	2271.47824	-0.4272710	0.0481656	-8.870869
ENSMUSG00000022884.17 Eif4a2	11561.93282	-0.4719840	0.0534312	-8.833493
ENSMUSG00000001891.17 Ugp2	1364.63153	-0.5993403	0.0678779	-8.829678
ENSMUSG00000084838.4 Gm10241	29.45661	8.3716581	0.9510624	8.802428
ENSMUSG00000026775.10 Yme1l1	1252.73281	-0.4786876	0.0544396	-8.792996
ENSMUSG00000020650.16 Bcap29	663.57029	-0.7001171	0.0803572	-8.712561
ENSMUSG00000042303.20 Sgsm3	1281.38919	0.3735305	0.0430203	8.682650
ENSMUSG00000071533.12 Pcnp	2553.37357	-0.4279186	0.0500775	-8.545132
ENSMUSG00000020238.15 Ncln	1446.32940	0.4261658	0.0499872	8.525490
ENSMUSG00000014905.5 Dnajb9	663.84669	-0.6632246	0.0778235	-8.522161
ENSMUSG00000057963.10 Itpk1	1451.17452	0.4522885	0.0532812	8.488701

mymds(de=dge3_o,name="Cont3: Effect of OVA treatment in hippocampus",myss=ss3_o,mx=xx3_o)

maplot(dge3_o,"Cont3: Effect of OVA treatment in hippocampus")

make_volcano(dge3_o,"Cont3: Effect of OVA treatment in hippocampus")

make_heatmap(de=dge3_o,name="Cont3: Effect of OVA treatment in hippocampus",myss=ss3_o,mx=xx3_o,n=50)

make_heatmap2(de=dge3_o,name="Cont3: Effect of OVA treatment in hippocampus",myss=ss3_o,mx=xx3_o,n=50)

dge3_o[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 3: Effect of treatment in hippocampus") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 3: Effect of treatment in hippocampus
	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000008668.16 Rps18	286.52172	-1.2451711	0.1674818	-7.434667	0.0e+00	0.0000000
ENSMUSG00000069008.4 Gm5537	14.35881	21.8370341	2.9884331	7.307185	0.0e+00	NA
ENSMUSG00000103034.2 Gm8797	107.80783	-1.1733731	0.1762296	-6.658206	0.0e+00	NA
ENSMUSG00000100862.2 Gm10925	1945.80382	5.8231351	0.9533584	6.108023	0.0e+00	0.0000059
ENSMUSG00000040447.16 Spns2	1527.35914	0.4988942	0.0886928	5.624970	0.0e+00	0.0000723
ENSMUSG00000081603.2 Gm14681	113.88393	-0.8342655	0.1522504	-5.479560	0.0e+00	0.0001247
ENSMUSG00000031858.17 Mau2	1851.99608	0.4053541	0.0759055	5.340248	1.0e-07	0.0002171
ENSMUSG00000084838.4 Gm10241	17.62940	7.6366750	1.4333637	5.327800	1.0e-07	NA
ENSMUSG00000021270.14 Hsp90aa1	19731.71056	-0.4150097	0.0812580	-5.107311	3.0e-07	0.0006369
ENSMUSG00000024883.9 Rin1	1730.36849	0.3961802	0.0785513	5.043585	5.0e-07	0.0007633
ENSMUSG00000001467.6 Cyp51	914.60992	-0.3266949	0.0652770	-5.004751	6.0e-07	0.0008176
ENSMUSG00000102627.2 Snurf	22.95625	8.0174386	1.6395147	4.890129	1.0e-06	NA
ENSMUSG00000044285.9 Ubb-ps	2010.40856	-0.9247134	0.1896289	-4.876438	1.1e-06	0.0014035
ENSMUSG00000107499.2 Ccdc142	44.66779	1.0467235	0.2150830	4.866603	1.1e-06	NA
ENSMUSG00000060862.11 Zbtb40	211.85127	0.6799751	0.1415661	4.803233	1.6e-06	0.0018257
ENSMUSG00000053550.14 Shisa7	2819.07578	0.6097051	0.1286592	4.738915	2.1e-06	0.0021791
ENSMUSG00000048154.18 Kmt2d	2495.07018	0.7784121	0.1651030	4.714707	2.4e-06	0.0021791
ENSMUSG00000032855.7 Pkd1	2513.84342	0.7124951	0.1511267	4.714553	2.4e-06	0.0021791
ENSMUSG00000078881.11 Gm14434	265.26689	-1.1328384	0.2439904	-4.642963	3.4e-06	0.0028688
ENSMUSG00000094437.3 Gm9830	43.94751	1.7576897	0.3813574	4.609036	4.0e-06	NA

Results for PAG

mymds(de=dge4_l,name="Cont4: Effect of LPS treatment in PAG",myss=ss4_l,mx=xx4_l)

maplot(dge4_l,"Cont4: Effect of LPS treatment in PAG")

make_volcano(dge4_l,"Cont4: Effect of LPS treatment in PAG")

make_heatmap(de=dge4_l,name="Cont4: Effect of LPS treatment in PAG",myss=ss4_l,mx=xx4_l,n=50)

make_heatmap2(de=dge4_l,name="Cont3: Effect of LPS treatment in PAG",myss=ss4_l,mx=xx4_l,n=50)

dge4_l[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 4: Effect of LPS treatment in PAG") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 4: Effect of LPS treatment in PAG
	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000101111.2 Gm28437	1489.59991	-28.4058069	2.9866471	-9.510935	0.0000000	0.0000000
ENSMUSG00000034459.9 Ifit1	55.64874	-1.1808759	0.1686328	-7.002649	0.0000000	0.0000000
ENSMUSG00000079363.8 Gbp4	95.91494	-1.2292487	0.1854179	-6.629613	0.0000000	0.0000002
ENSMUSG00000025498.16 Irf7	84.22641	-1.0483087	0.1595935	-6.568618	0.0000000	0.0000002
ENSMUSG00000078922.10 Tgtp1	21.82329	-2.6093058	0.4065209	-6.418626	0.0000000	0.0000005
ENSMUSG00000102070.2 Gm28661	10196.94116	-14.7763448	2.9540986	-5.001981	0.0000006	0.0016656
ENSMUSG00000028820.14 Sfpq	2198.16793	0.2648449	0.0535496	4.945791	0.0000008	0.0019080
ENSMUSG00000036915.18 Kirrel2	82.56463	0.9736571	0.2004918	4.856344	0.0000012	0.0026324
ENSMUSG00000026222.17 Sp100	58.86066	-0.8315047	0.1797559	-4.625743	0.0000037	0.0073042
ENSMUSG00000075602.11 Ly6a	107.28664	-0.8323082	0.1822514	-4.566814	0.0000050	0.0087213
ENSMUSG00000024079.5 Eif2ak2	182.29278	-0.3940654	0.0893055	-4.412555	0.0000102	0.0155853
ENSMUSG00000078921.4 Tgtp2	20.19923	-1.2810520	0.2908727	-4.404167	0.0000106	0.0155853
ENSMUSG00000052776.11 Oas1a	19.16469	-1.2619796	0.2896101	-4.357512	0.0000132	0.0178219
ENSMUSG00000028268.15 Gbp3	138.45727	-0.5635390	0.1327639	-4.244671	0.0000219	0.0275394
ENSMUSG00000054072.13 Iigp1	23.29794	-1.4762784	0.3644390	-4.050824	0.0000510	0.0599248
ENSMUSG00000110862.2 Gm35835	17.20473	-1.2465344	0.3123628	-3.990663	0.0000659	0.0685953
ENSMUSG00000027875.13 Hmgcs2	123.00353	-0.7999616	0.2005166	-3.989504	0.0000662	0.0685953
ENSMUSG00000002289.17 Angptl4	67.16713	-0.6692079	0.1702196	-3.931438	0.0000844	0.0826191
ENSMUSG00000015085.9 Entpd2	359.15806	0.3576188	0.0921084	3.882584	0.0001034	0.0954024
ENSMUSG00000030107.11 Usp18	27.77358	-0.9172425	0.2372305	-3.866461	0.0001104	0.0954024

mymds(de=dge4_o,name="Cont4: Effect of OVA treatment in PAG",myss=ss4_o,mx=xx4_o)

maplot(dge4_o,"Cont4: Effect of OVA treatment in PAG")

make_volcano(dge4_o,"Cont4: Effect of OVA treatment in PAG")

make_heatmap(de=dge4_o,name="Cont4: Effect of OVA treatment in PAG",myss=ss4_o,mx=xx4_o,n=50)

make_heatmap2(de=dge4_o,name="Cont4: Effect of OVA treatment in PAG",myss=ss4_o,mx=xx4_o,n=50)

dge4_o[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 4: Effect of treatment in PAG") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 4: Effect of treatment in PAG
	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000101111.2 Gm28437	1504.44501	-28.8708740	2.9866471	-9.666651	0.0000000	0.0000000
ENSMUSG00000036915.18 Kirrel2	88.57030	1.1004532	0.2077438	5.297165	0.0000001	0.0010387
ENSMUSG00000106604.3 Gm42535	45.90967	-1.3651174	0.2757149	-4.951192	0.0000007	0.0043425
ENSMUSG00000102070.2 Gm28661	10298.72248	-14.1414026	2.9291777	-4.827772	0.0000014	0.0060964
ENSMUSG00000085328.3 Gm17131	231.90596	-0.6897420	0.1730516	-3.985758	0.0000673	0.2376070
ENSMUSG00000118642.2 AY036118	813.47892	-0.9127923	0.2465059	-3.702922	0.0002131	0.6273842
ENSMUSG00000102289.2 Gm31258	14.36741	-1.5814526	0.4320332	-3.660488	0.0002517	0.6351635
ENSMUSG00000021906.15 Oxnad1	298.24353	0.2627823	0.0731645	3.591666	0.0003286	0.7254015
ENSMUSG00000028820.14 Sfpq	2221.70651	0.2671342	0.0763080	3.500736	0.0004640	0.9105264
ENSMUSG00000094672.3 Vmn2r25	14.38564	-1.6452102	0.4753780	-3.460847	0.0005385	0.9510627
ENSMUSG00000054252.19 Fgfr3	2567.78162	0.3560764	0.1040765	3.421295	0.0006232	0.9754406
ENSMUSG00000024998.18 Plce1	525.21275	0.3861609	0.1134704	3.403185	0.0006661	0.9754406
ENSMUSG00000013089.16 Etv5	501.24250	0.2879667	0.0851978	3.379980	0.0007249	0.9754406
ENSMUSG00000047216.9 Cdh19	172.79548	0.5469923	0.1627061	3.361844	0.0007742	0.9754406
ENSMUSG00000115632.2 Gm17922	12.17839	-2.0478089	0.6150693	-3.329396	0.0008703	0.9754406
ENSMUSG00000067220.13 Cnga1	22.07693	-1.3574772	0.4082491	-3.325120	0.0008838	0.9754406
ENSMUSG00000046678.7 Olfr981	15.67008	-1.2921922	0.3905994	-3.308229	0.0009389	0.9754406
ENSMUSG00000007682.8 Dio2	395.46383	0.4716641	0.1435644	3.285383	0.0010184	0.9904947
ENSMUSG00000108827.4 Olfr1310	28.09654	-1.3088894	0.4003901	-3.269035	0.0010791	0.9904947
ENSMUSG00000002831.15 Plin4	203.85756	1.4385797	0.4416960	3.256946	0.0011262	0.9904947

Results for NI

mymds(de=dge5_l,name="Cont5: Effect of LPS treatment in NI",myss=ss5_l,mx=xx5_l)

maplot(dge5_l,"Cont5: Effect of LPS treatment in NI")

make_volcano(dge5_l,"Cont5: Effect of LPS treatment in NI")

make_heatmap(de=dge5_l,name="Cont5: Effect of LPS treatment in NI",myss=ss5_l,mx=xx5_l,n=50)

make_heatmap2(de=dge5_l,name="Cont3: Effect of LPS treatment in NI",myss=ss5_l,mx=xx5_l,n=50)

dge5_l[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 5: Effect of LPS treatment in NI") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 5: Effect of LPS treatment in NI
	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000004655.6 Aqp1	64.33666	-2.5245691	0.4256671	-5.930853	0.0e+00	0.0000404
ENSMUSG00000002831.15 Plin4	250.84384	1.5543801	0.2773464	5.604472	0.0e+00	0.0001368
ENSMUSG00000026822.15 Lcn2	16.22880	-2.6059016	0.4734932	-5.503567	0.0e+00	NA
ENSMUSG00000028919.12 Arhgef19	395.89057	0.4100393	0.0746763	5.490892	0.0e+00	0.0001368
ENSMUSG00000022425.17 Enpp2	8280.82063	-1.1048917	0.2013475	-5.487486	0.0e+00	0.0001368
ENSMUSG00000068699.13 Flnc	125.76677	0.8058438	0.1494388	5.392467	1.0e-07	0.0001573
ENSMUSG00000053646.14 Plxnb1	2389.38786	0.3753821	0.0697524	5.381640	1.0e-07	0.0001573
ENSMUSG00000071267.12 Zfp942	226.23109	-0.5520131	0.1029389	-5.362533	1.0e-07	0.0001573
ENSMUSG00000036890.14 Gtdc1	411.57310	-0.4107206	0.0781476	-5.255700	1.0e-07	0.0002474
ENSMUSG00000078922.10 Tgtp1	37.06696	-2.2129060	0.4225152	-5.237459	2.0e-07	NA
ENSMUSG00000039959.14 Hip1	856.67943	0.5239050	0.1013262	5.170477	2.0e-07	0.0003482
ENSMUSG00000068323.13 Slc4a5	92.61322	-1.8988675	0.3702782	-5.128219	3.0e-07	0.0003771
ENSMUSG00000114133.2 Gm20075	468.65132	-0.3716223	0.0726130	-5.117845	3.0e-07	0.0003771
ENSMUSG00000021696.10 Elovl7	524.66657	-0.4441318	0.0873935	-5.081979	4.0e-07	0.0004178
ENSMUSG00000027875.13 Hmgcs2	132.26542	-1.1687659	0.2334476	-5.006545	6.0e-07	0.0005721
ENSMUSG00000060862.11 Zbtb40	245.30142	0.6788979	0.1360226	4.991069	6.0e-07	0.0005756
ENSMUSG00000113337.2 Gm19220	2324.88228	1.0465161	0.2131860	4.908936	9.0e-07	0.0008193
ENSMUSG00000045039.10 Megf8	3979.39321	0.3979746	0.0822716	4.837328	1.3e-06	0.0011038
ENSMUSG00000069833.14 Ahnak	814.46670	0.6652302	0.1391700	4.779984	1.8e-06	0.0013840
ENSMUSG00000015652.10 Steap1	18.65542	-2.4798774	0.5203028	-4.766219	1.9e-06	NA

mymds(de=dge5_o,name="Cont5: Effect of OVA treatment in NI",myss=ss5_o,mx=xx5_o)

maplot(dge5_o,"Cont4: Effect of OVA treatment in NI")

make_volcano(dge5_o,"Cont5: Effect of OVA treatment in NI")

make_heatmap(de=dge5_o,name="Cont5: Effect of OVA treatment in NI",myss=ss5_o,mx=xx5_o,n=50)

make_heatmap2(de=dge5_o,name="Cont5: Effect of OVA treatment in NI",myss=ss5_o,mx=xx5_o,n=50)

dge5_o[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 5: Effect of treatment in NI") %>%
  kable_paper("hover", full_width = F)

Top gene expression differences for contrast 5: Effect of treatment in NI
	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000056260.16 Lrif1	211.98587	-0.7336887	0.1100875	-6.664595	0.00e+00	0.0000003
ENSMUSG00000002831.15 Plin4	258.40488	1.6709594	0.2980839	5.605668	0.00e+00	0.0001173
ENSMUSG00000055322.17 Tns1	1466.90156	0.3981360	0.0721780	5.516031	0.00e+00	0.0001307
ENSMUSG00000015944.10 Castor2	3213.72353	0.3932529	0.0758616	5.183820	2.00e-07	0.0006144
ENSMUSG00000024175.3 Tekt4	67.10906	0.9393058	0.1822044	5.155233	3.00e-07	NA
ENSMUSG00000020116.8 Pno1	507.57144	-0.3892492	0.0798403	-4.875347	1.10e-06	0.0024561
ENSMUSG00000000194.14 Gpr107	1069.00316	0.3266125	0.0685118	4.767246	1.90e-06	0.0035192
ENSMUSG00000055782.10 Abcd2	665.20587	-0.4696031	0.0996124	-4.714301	2.40e-06	0.0035595
ENSMUSG00000117284.2 Gm7072	208.57068	-0.5445045	0.1156892	-4.706613	2.50e-06	0.0035595
ENSMUSG00000074505.7 Fat3	548.94128	0.5489851	0.1177974	4.660416	3.20e-06	0.0039643
ENSMUSG00000060862.11 Zbtb40	233.67044	0.5799773	0.1251034	4.635981	3.60e-06	0.0040164
ENSMUSG00000069833.14 Ahnak	857.16074	0.8065828	0.1763222	4.574483	4.80e-06	0.0046838
ENSMUSG00000067038.7 Rps12-ps3	447.47971	-1.9169288	0.4198274	-4.565993	5.00e-06	0.0046838
ENSMUSG00000027708.15 Dcun1d1	788.59922	-0.3874162	0.0869258	-4.456861	8.30e-06	0.0072331
ENSMUSG00000045039.10 Megf8	3981.85021	0.4012817	0.0919057	4.366235	1.26e-05	0.0101992
ENSMUSG00000043733.15 Ptpn11	4319.80686	0.3173738	0.0730162	4.346623	1.38e-05	0.0101992
ENSMUSG00000019232.15 Etnppl	1154.63560	0.7040907	0.1623392	4.337157	1.44e-05	0.0101992
ENSMUSG00000068798.11 Rap1a	1471.01059	-0.4872519	0.1137471	-4.283643	1.84e-05	0.0115374
ENSMUSG00000027206.14 Cops2	2031.84090	-0.3788991	0.0885072	-4.280998	1.86e-05	0.0115374
ENSMUSG00000025903.15 Lypla1	715.48002	-0.5016851	0.1174408	-4.271815	1.94e-05	0.0115374

Euler diagrams

Below we’re comparing the effect of the two treatments on the number of genes found, using an Euler diagram.

par(mar=c(2,2,2,2))

v0 <- list("hyp LPS up"=d1up_l, 
  "hyp LPS dn"=d1dn_l,
  "hyp OVA up"=d1up_o,
  "hyp OVA dn"=d1dn_o)
plot(euler(v0),quantities = TRUE, edges = "gray", main="effect of treatments on hyp")

v0 <- list("amy LPS up"=d2up_l,
  "amy LPS dn"=d2dn_l,
  "amy OVA up"=d2up_o,
  "amy OVA dn"=d2dn_o)
plot(euler(v0),quantities = TRUE, edges = "gray", main="effect of treatments on amy")

v0 <- list("hip LPS up"=d3up_l,
  "hip LPS dn"=d3dn_l,
  "hip OVA up"=d3up_o,
  "hip OVA dn"=d3dn_o)
plot(euler(v0),quantities = TRUE, edges = "gray", main="effect of treatments on hip")

v0 <- list("pag LPS up"=d4up_l,
  "pag LPS dn"=d4dn_l,
  "pag OVA up"=d4up_o,
  "pag OVA dn"=d4dn_o)
plot(euler(v0),quantities = TRUE, edges = "gray", main="effect of treatments on pag")

v0 <- list("NI LPS up"=d5up_l,
  "NI LPS dn"=d5dn_l,
  "NI OVA up"=d5up_o,
  "NI OVA dn"=d5dn_o)
plot(euler(v0),quantities = TRUE, edges = "gray", main="effect of treatments on ni")

Mitch pathway analysis in hyp

Gene sets are obtained from Reactome.

Firstly need to conduct mitch enrichment analysis for each contrast separately.

if  (!file.exists("mouse_msigdb_reactome_2022-02-16.gmt")) {
  download.file("http://ziemann-lab.net/public/msigdb_mouse/mouse_msigdb_reactome_2022-02-16.gmt",
    destfile="mouse_msigdb_reactome_2022-02-16.gmt")
}
genesets <- gmt_import("mouse_msigdb_reactome_2022-02-16.gmt")
names(genesets) <- gsub("REACTOME_","",names(genesets))
names(genesets) <- gsub("_"," ",names(genesets))

# gene table
gt <- as.data.frame(rownames(xx))
gt$gene <- sapply(strsplit(gt[,1]," "),"[[",2)

Mitch for LPS in hyp

m1_l <- mitch_import(dge1_l, DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 17095

## Note: no. genes in output = 17078

## Note: estimated proportion of input genes in output = 0.999

mres1_l <- mitch_calc(m1_l, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres1_l$enrichment_result,20) %>% 
  kbl(caption = "Top gene pathway differences in hypothalamus for LPS") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in hypothalamus for LPS
	set	setSize	pANOVA	s.dist	p.adjustANOVA
726	POST CHAPERONIN TUBULIN FOLDING PATHWAY	17	0.0000074	0.6275717	0.0005720
1138	TRANSPORT OF CONNEXONS TO THE PLASMA MEMBRANE	13	0.0002277	0.5904077	0.0059649
337	FORMATION OF TUBULIN FOLDING INTERMEDIATES BY CCT TRIC	19	0.0000123	0.5793238	0.0007249
143	CGMP EFFECTS	15	0.0001524	-0.5647268	0.0047278
785	REDUCTION OF CYTOSOLIC CA LEVELS	11	0.0013336	-0.5586911	0.0212474
629	NITRIC OXIDE STIMULATES GUANYLATE CYCLASE	20	0.0000599	-0.5183198	0.0022099
910	SEALING OF THE NUCLEAR ENVELOPE NE BY ESCRT III	23	0.0000223	0.5107134	0.0010992
679	P75NTR REGULATES AXONOGENESIS	10	0.0068797	-0.4935552	0.0659004
66	APOPTOSIS INDUCED DNA FRAGMENTATION	10	0.0091035	-0.4763065	0.0800971
1173	WNT5A DEPENDENT INTERNALIZATION OF FZD2 FZD5 AND ROR2	13	0.0033723	0.4695846	0.0401608
161	CLEC7A DECTIN 1 INDUCES NFAT ACTIVATION	11	0.0108084	-0.4438283	0.0912425
368	GAP JUNCTION ASSEMBLY	21	0.0005678	0.4344453	0.0111580
498	IRE1ALPHA ACTIVATES CHAPERONES	50	0.0000004	0.4142824	0.0000794
388	GLYCOGEN STORAGE DISEASES	12	0.0142457	0.4086390	0.1072601
1170	VOLTAGE GATED POTASSIUM CHANNELS	39	0.0000143	-0.4015720	0.0008015
112	CARBOXYTERMINAL POST TRANSLATIONAL MODIFICATIONS OF TUBULIN	36	0.0000416	0.3947013	0.0018179
43	AGGREPHAGY	35	0.0000581	0.3926790	0.0022099
1072	THE ROLE OF GTSE1 IN G2 M PROGRESSION AFTER G2 CHECKPOINT	66	0.0000001	0.3869194	0.0000214
201	CYTOSOLIC IRON SULFUR CLUSTER ASSEMBLY	13	0.0163463	0.3846244	0.1160336
1006	SODIUM CALCIUM EXCHANGERS	10	0.0353060	-0.3844036	0.1963481

m1_l_top <- subset(mres1_l$enrichment_result,p.adjustANOVA<0.05)
m1_l_up <- subset(m1_l_top,s.dist>0)$set
m1_l_dn <- subset(m1_l_top,s.dist<0)$set
#mitch_report(mres1_l,outfile="mitch1_l.html",overwrite=TRUE)
write.table(mres1_l$enrichment_result,file="mitch1_l.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m1_l_top_up <- head(subset(m1_l_top,s.dist>0),10)[,"s.dist"]
names(m1_l_top_up) <- head(subset(m1_l_top,s.dist>0),10)[,"set"]
m1_l_top_dn <- head(subset(m1_l_top,s.dist<0),10)[,"s.dist"]
names(m1_l_top_dn) <- head(subset(m1_l_top,s.dist<0),10)[,"set"]
m1_l_top_updn <- c(m1_l_top_up,m1_l_top_dn)
m1_l_top_updn <- m1_l_top_updn[order(m1_l_top_updn)]

par(mar=c(5,25,3,1))
barplot(m1_l_top_updn,horiz=TRUE,las=1,col="darkgray",
  xlab="Enrichment score",cex.names=0.8,xlim=c(-1,1),
  main="Effect of LPS in hyp")
grid()

Mitch for OVA in hyp

m1_o <- mitch_import(dge1_o, DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 17114

## Note: no. genes in output = 17098

## Note: estimated proportion of input genes in output = 0.999

mres1_o <- mitch_calc(m1_o, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres1_o$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in hypothalamus for OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in hypothalamus for OVA
	set	setSize	pANOVA	s.dist	p.adjustANOVA
304	EUKARYOTIC TRANSLATION ELONGATION	87	0.0000000	0.7529442	0.0000000
1011	SRP DEPENDENT COTRANSLATIONAL PROTEIN TARGETING TO MEMBRANE	106	0.0000000	0.6728738	0.0000000
838	RESPONSE OF EIF2AK4 GCN2 TO AMINO ACID DEFICIENCY	95	0.0000000	0.6489214	0.0000000
305	EUKARYOTIC TRANSLATION INITIATION	114	0.0000000	0.6266445	0.0000000
790	REGULATION OF COMMISSURAL AXON PATHFINDING BY SLIT AND ROBO	10	0.0006068	-0.6260885	0.0066803
1175	YAP1 AND WWTR1 TAZ STIMULATED GENE EXPRESSION	11	0.0004339	0.6126348	0.0050606
28	ACTIVATION OF THE MRNA UPON BINDING OF THE CAP BINDING COMPLEX AND EIFS AND SUBSEQUENT BINDING TO 43S	59	0.0000000	0.5829488	0.0000000
912	SELENOAMINO ACID METABOLISM	109	0.0000000	0.5767461	0.0000000
27	ACTIVATION OF THE AP 1 FAMILY OF TRANSCRIPTION FACTORS	10	0.0024530	-0.5531601	0.0209392
328	FORMATION OF ATP BY CHEMIOSMOTIC COUPLING	18	0.0000765	0.5384075	0.0013456
631	NONSENSE MEDIATED DECAY NMD	109	0.0000000	0.5242653	0.0000000
171	COMPLEX I BIOGENESIS	56	0.0000000	0.5180541	0.0000000
157	CLASS C 3 METABOTROPIC GLUTAMATE PHEROMONE RECEPTORS	12	0.0019791	-0.5157049	0.0175292
923	SEROTONIN NEUROTRANSMITTER RELEASE CYCLE	18	0.0001961	-0.5069737	0.0028334
8	ACETYLCHOLINE BINDING AND DOWNSTREAM EVENTS	10	0.0061207	-0.5006086	0.0450635
419	HIGHLY CALCIUM PERMEABLE POSTSYNAPTIC NICOTINIC ACETYLCHOLINE RECEPTORS	10	0.0061207	-0.5006086	0.0450635
283	ENDOSOMAL VACUOLAR PATHWAY	11	0.0048523	0.4904579	0.0381065
1049	SYNTHESIS OF PROSTAGLANDINS PG AND THROMBOXANES TX	13	0.0031670	0.4726998	0.0265946
791	REGULATION OF EXPRESSION OF SLITS AND ROBOS	160	0.0000000	0.4702680	0.0000000
836	RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CHEMIOSMOTIC COUPLING AND HEAT PRODUCTION BY UNCOUPLING PROTEINS	125	0.0000000	0.4602665	0.0000000

m1_o_top <- subset(mres1_o$enrichment_result,p.adjustANOVA<0.05)
m1_o_up <- subset(m1_o_top,s.dist>0)$set
m1_o_dn <- subset(m1_o_top,s.dist<0)$set
#mitch_report(mres1_o,outfile="mitch1_o.html",overwrite=TRUE)
write.table(mres1_o$enrichment_result,file="mitch1_o.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m1_o_top_up <- head(subset(m1_o_top,s.dist>0),10)[,"s.dist"]
names(m1_o_top_up) <- head(subset(m1_o_top,s.dist>0),10)[,"set"]
m1_o_top_dn <- head(subset(m1_o_top,s.dist<0),10)[,"s.dist"]
names(m1_o_top_dn) <- head(subset(m1_o_top,s.dist<0),10)[,"set"]
m1_o_top_updn <- c(m1_o_top_up,m1_o_top_dn)
m1_o_top_updn <- m1_o_top_updn[order(m1_o_top_updn)]

par(mar=c(5,25,3,1))
barplot(m1_o_top_updn,horiz=TRUE,las=1,col="darkgray",
  xlab="Enrichment score",cex.names=0.8,xlim=c(-1,1),
  main="Effect of LPS in hyp")
grid()

Two dimensional mitch analysis in hyp.

ml <- list("LPS"=dge1_l,"OVA"=dge1_o)
m1 <- mitch_import(ml, DEtype="deseq2",geneTable=gt)

## Note: Mean no. genes in input = 17104.5

## Note: no. genes in output = 16956

## Note: estimated proportion of input genes in output = 0.991

head(m1)

##                      LPS         OVA
## 0610009B22Rik -0.3304181  0.09722066
## 0610009E02Rik  0.7596333  1.26326253
## 0610009L18Rik -1.2521412 -0.66374969
## 0610010K14Rik  0.4193627 -0.45368981
## 0610012G03Rik  1.4462155  0.12236243
## 0610030E20Rik -0.2426628 -0.04562294

mres1 <- mitch_calc(m1, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(mres1$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in hypothalamus for LPS and OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in hypothalamus for LPS and OVA
	set	setSize	pMANOVA	s.LPS	s.OVA	p.LPS	p.OVA	s.dist	SD	p.adjustMANOVA
304	EUKARYOTIC TRANSLATION ELONGATION	87	0.0000000	0.3003667	0.7570365	0.0000013	0.0000000	0.8144473	0.3229143	0.0000000
1008	SRP DEPENDENT COTRANSLATIONAL PROTEIN TARGETING TO MEMBRANE	106	0.0000000	0.2540227	0.6770203	0.0000063	0.0000000	0.7231072	0.2991045	0.0000000
837	RESPONSE OF EIF2AK4 GCN2 TO AMINO ACID DEFICIENCY	95	0.0000000	0.2222176	0.6529768	0.0001830	0.0000000	0.6897531	0.3045928	0.0000000
789	REGULATION OF COMMISSURAL AXON PATHFINDING BY SLIT AND ROBO	10	0.0029173	-0.2597191	-0.6239112	0.1550094	0.0006340	0.6758101	0.2575227	0.0211031
305	EUKARYOTIC TRANSLATION INITIATION	114	0.0000000	0.2043357	0.6307946	0.0001655	0.0000000	0.6630647	0.3015520	0.0000000
783	REDUCTION OF CYTOSOLIC CA LEVELS	11	0.0047499	-0.5569302	-0.3356045	0.0013813	0.0539498	0.6502320	0.1565009	0.0308348
1172	YAP1 AND WWTR1 TAZ STIMULATED GENE EXPRESSION	11	0.0015182	0.1440328	0.6163095	0.4081891	0.0004006	0.6329161	0.3339500	0.0127425
724	POST CHAPERONIN TUBULIN FOLDING PATHWAY	17	0.0000053	0.6288482	-0.0037817	0.0000071	0.9784665	0.6288595	0.4473369	0.0001307
910	SELENOAMINO ACID METABOLISM	109	0.0000000	0.2406957	0.5806511	0.0000143	0.0000000	0.6285620	0.2403848	0.0000000
337	FORMATION OF TUBULIN FOLDING INTERMEDIATES BY CCT TRIC	19	0.0000643	0.5807497	0.1992927	0.0000117	0.1326527	0.6139933	0.2697308	0.0010071
28	ACTIVATION OF THE MRNA UPON BINDING OF THE CAP BINDING COMPLEX AND EIFS AND SUBSEQUENT BINDING TO 43S	59	0.0000000	0.1619654	0.5869250	0.0314726	0.0000000	0.6088626	0.3004918	0.0000000
143	CGMP EFFECTS	15	0.0008040	-0.5628515	-0.2196997	0.0001603	0.1407302	0.6042100	0.2426449	0.0080064
1135	TRANSPORT OF CONNEXONS TO THE PLASMA MEMBRANE	13	0.0007836	0.5915536	0.1216795	0.0002214	0.4475253	0.6039384	0.3322512	0.0078692
677	P75NTR REGULATES AXONOGENESIS	10	0.0197324	-0.4916322	-0.3297415	0.0071009	0.0709981	0.5919728	0.1144740	0.0916424
27	ACTIVATION OF THE AP 1 FAMILY OF TRANSCRIPTION FACTORS	10	0.0104603	-0.2020772	-0.5509737	0.2685386	0.0025521	0.5868622	0.2467071	0.0582503
618	NEUROTOXICITY OF CLOSTRIDIUM TOXINS	10	0.0247953	-0.3756757	-0.4495692	0.0396833	0.0138279	0.5858709	0.0522506	0.1067196
410	HDR THROUGH MMEJ ALT NHEJ	10	0.0287449	-0.3548094	-0.4485542	0.0520448	0.0140436	0.5719184	0.0662876	0.1189268
283	ENDOSOMAL VACUOLAR PATHWAY	11	0.0001934	-0.2756673	0.4953298	0.1134178	0.0044459	0.5668722	0.5451773	0.0024170
157	CLASS C 3 METABOTROPIC GLUTAMATE PHEROMONE RECEPTORS	12	0.0086541	-0.2343602	-0.5130823	0.1598426	0.0020867	0.5640729	0.1970863	0.0496027
631	NONSENSE MEDIATED DECAY NMD	109	0.0000000	0.1971761	0.5281756	0.0003787	0.0000000	0.5637800	0.2340520	0.0000000

m1_top <- subset(mres1$enrichment_result,p.adjustMANOVA<0.05)
m1_up <- subset(m1_top,s.dist>0)$set
m1_dn <- subset(m1_top,s.dist<0)$set
#mitch_report(mres1,outfile="mitch1.html",overwrite=TRUE)
write.table(mres1$enrichment_result,file="mitch1.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m1_top_up <- head(subset(m1_top,s.dist>0),10)[,"s.dist"]
names(m1_top_up) <- head(subset(m1_top,s.dist>0),10)[,"set"]
m1_top_dn <- head(subset(m1_top,s.dist<0),10)[,"s.dist"]
names(m1_top_dn) <- head(subset(m1_top,s.dist<0),10)[,"set"]
m1_top_updn <- c(m1_top_up,m1_top_dn)
m1_top_updn <- m1_top_updn[order(m1_top_updn)]

# heatmap
mr1 <- head(mres1$enrichment_result,30)
mr1 <- as.matrix(mr1[,4:5])
rownames(mr1) <- head(mres1$enrichment_result,30)$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mr1,trace="none",col=colfunc(25),scale="none", margins = c(10,28), cexRow=0.7,
  main="Top ranked pathways in hypothalamus" , cexCol=0.9 )

Mitch pathway analysis in amygdala

Mitch for LPS in amy

m2_l <- mitch_import(dge2_l, DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 16818

## Note: no. genes in output = 16805

## Note: estimated proportion of input genes in output = 0.999

mres2_l <- mitch_calc(m2_l, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres2_l$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in amy for LPS") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in amy for LPS
	set	setSize	pANOVA	s.dist	p.adjustANOVA
534	MET ACTIVATES RAP1 AND RAC1	11	0.0010392	-0.5710589	0.0295794
1039	SYNTHESIS OF PIPS AT THE LATE ENDOSOME MEMBRANE	11	0.0021453	-0.5344441	0.0510929
24	ACTIVATION OF RAC1	13	0.0017643	-0.5009253	0.0447585
443	INITIAL TRIGGERING OF COMPLEMENT	10	0.0100369	0.4701756	0.1541198
199	CYTOSOLIC IRON SULFUR CLUSTER ASSEMBLY	13	0.0063266	0.4373603	0.1153611
301	EUKARYOTIC TRANSLATION ELONGATION	87	0.0000000	0.4368724	0.0000000
185	CREB1 PHOSPHORYLATION THROUGH THE ACTIVATION OF ADENYLATE CYCLASE	11	0.0121751	-0.4365520	0.1711852
621	NF KB IS ACTIVATED AND SIGNALS SURVIVAL	13	0.0068114	0.4334482	0.1222900
20	ACTIVATION OF IRF3 IRF7 MEDIATED BY TBK1 IKK EPSILON	16	0.0027795	0.4319123	0.0636015
909	SEMA4D INDUCED CELL MIGRATION AND GROWTH CONE COLLAPSE	20	0.0009094	0.4284599	0.0265313
901	SCAVENGING BY CLASS A RECEPTORS	14	0.0075549	0.4123553	0.1306844
917	SEROTONIN RECEPTORS	10	0.0291904	-0.3982971	0.2958961
164	COLLAGEN CHAIN TRIMERIZATION	40	0.0000146	0.3960722	0.0008119
1037	SYNTHESIS OF PIPS AT THE EARLY ENDOSOME MEMBRANE	16	0.0061523	-0.3955938	0.1139641
798	REGULATION OF INNATE IMMUNE RESPONSES TO CYTOSOLIC DNA	13	0.0188749	0.3761223	0.2265723
833	RESPONSE OF EIF2AK4 GCN2 TO AMINO ACID DEFICIENCY	94	0.0000000	0.3738816	0.0000001
45	ALPHA LINOLENIC OMEGA3 AND LINOLEIC OMEGA6 ACID METABOLISM	12	0.0277431	0.3669485	0.2916776
163	COLLAGEN BIOSYNTHESIS AND MODIFYING ENZYMES	63	0.0000010	0.3568954	0.0000704
629	NONSENSE MEDIATED DECAY NMD	109	0.0000000	0.3466105	0.0000001
42	AGGREPHAGY	35	0.0005204	0.3389590	0.0164130

m2_l_top <- subset(mres2_l$enrichment_result,p.adjustANOVA<0.05)
m2_l_up <- subset(m2_l_top,s.dist>0)$set
m2_l_dn <- subset(m2_l_top,s.dist<0)$set
#mitch_report(mres2_l,outfile="mitch2_l.html",overwrite=TRUE)
write.table(mres2_l$enrichment_result,file="mitch2_l.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m2_l_top_up <- head(subset(m2_l_top,s.dist>0),10)[,"s.dist"]
names(m2_l_top_up) <- head(subset(m2_l_top,s.dist>0),10)[,"set"]
m2_l_top_dn <- head(subset(m2_l_top,s.dist<0),10)[,"s.dist"]
names(m2_l_top_dn) <- head(subset(m2_l_top,s.dist<0),10)[,"set"]
m2_l_top_updn <- c(m2_l_top_up,m2_l_top_dn)
m2_l_top_updn <- m2_l_top_updn[order(m2_l_top_updn)]

par(mar=c(5,25,3,1))
barplot(m2_l_top_updn,horiz=TRUE,las=1,col="darkgray",
  xlab="Enrichment score",cex.names=0.8,xlim=c(-1,1),
  main="Effect of LPS in amy")
grid()

Mitch for OVA in amy

m2_o <- mitch_import(dge2_o, DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 16624

## Note: no. genes in output = 16611

## Note: estimated proportion of input genes in output = 0.999

mres2_o <- mitch_calc(m2_o, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres2_o$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in amy for OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in amy for OVA
	set	setSize	pANOVA	s.dist	p.adjustANOVA
1096	TRAFFICKING AND PROCESSING OF ENDOSOMAL TLR	11	0.0001763	0.6530559	0.0016590
656	NUCLEOTIDE LIKE PURINERGIC RECEPTORS	13	0.0005347	0.5546451	0.0043939
39	ADP SIGNALLING THROUGH P2Y PURINOCEPTOR 12	19	0.0000407	0.5437560	0.0005103
72	ASPARTATE AND ASPARAGINE METABOLISM	10	0.0030231	0.5415336	0.0172096
170	COMPLEX I BIOGENESIS	56	0.0000000	0.5352440	0.0000000
1045	SYNTHESIS OF VERY LONG CHAIN FATTY ACYL COAS	20	0.0000351	0.5344042	0.0004546
1137	TRIGLYCERIDE CATABOLISM	14	0.0006681	0.5251638	0.0053182
704	PINK1 PRKN MEDIATED MITOPHAGY	22	0.0000244	0.5196982	0.0003353
344	FRS MEDIATED FGFR4 SIGNALING	12	0.0019164	0.5173002	0.0121544
921	SHC MEDIATED CASCADE FGFR4	10	0.0047384	0.5157762	0.0243601
302	EUKARYOTIC TRANSLATION ELONGATION	87	0.0000000	0.5144868	0.0000000
188	CRMPS IN SEMA3A SIGNALING	16	0.0003875	-0.5123682	0.0033247
831	RESPIRATORY ELECTRON TRANSPORT	102	0.0000000	0.5111484	0.0000000
153	CITRIC ACID CYCLE TCA CYCLE	22	0.0000354	0.5092805	0.0004546
700	PI 3K CASCADE FGFR4	10	0.0053424	0.5087163	0.0264908
832	RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CHEMIOSMOTIC COUPLING AND HEAT PRODUCTION BY UNCOUPLING PROTEINS	125	0.0000000	0.5039253	0.0000000
718	PLATELET SENSITIZATION BY LDL	15	0.0007660	0.5017755	0.0058008
343	FRS MEDIATED FGFR3 SIGNALING	14	0.0012050	0.4997547	0.0082239
950	SIGNALING BY FGFR3 FUSIONS IN CANCER	10	0.0062828	0.4990422	0.0302976
920	SHC MEDIATED CASCADE FGFR3	12	0.0029553	0.4955419	0.0169061

m2_o_top <- subset(mres2_o$enrichment_result,p.adjustANOVA<0.05)
m2_o_up <- subset(m2_o_top,s.dist>0)$set
m2_o_dn <- subset(m2_o_top,s.dist<0)$set
#mitch_report(mres2_o,outfile="mitch2_o.html",overwrite=TRUE)
write.table(mres2_o$enrichment_result,file="mitch2_o.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m2_o_top_up <- head(subset(m2_o_top,s.dist>0),10)[,"s.dist"]
names(m2_o_top_up) <- head(subset(m2_o_top,s.dist>0),10)[,"set"]
m2_o_top_dn <- head(subset(m2_o_top,s.dist<0),10)[,"s.dist"]
names(m2_o_top_dn) <- head(subset(m2_o_top,s.dist<0),10)[,"set"]
m2_o_top_updn <- c(m2_o_top_up,m2_o_top_dn)
m2_o_top_updn <- m2_o_top_updn[order(m2_o_top_updn)]

par(mar=c(5,25,3,1))
barplot(m2_o_top_updn,horiz=TRUE,las=1,col="darkgray",
  xlab="Enrichment score",cex.names=0.8,xlim=c(-1,1),
  main="Effect of LPS in amy")
grid()

Two dimensional mitch analysis in amy.

ml <- list("LPS"=dge2_l,"OVA"=dge2_o)
m2 <- mitch_import(ml, DEtype="deseq2",geneTable=gt)

## Note: Mean no. genes in input = 16721

## Note: no. genes in output = 16550

## Note: estimated proportion of input genes in output = 0.99

head(m2)

##                      LPS          OVA
## 0610009B22Rik  0.2684894  1.831236577
## 0610009E02Rik -0.7521235 -0.103702849
## 0610009L18Rik -1.1838923 -0.075637695
## 0610010K14Rik -0.0433030  0.417670579
## 0610012G03Rik -0.1930476  0.008709563
## 0610030E20Rik  1.0496464 -1.609523779

mres2 <- mitch_calc(m2, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(mres2$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in amy for LPS and OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in amy for LPS and OVA
	set	setSize	pMANOVA	s.LPS	s.OVA	p.LPS	p.OVA	s.dist	SD	p.adjustMANOVA
301	EUKARYOTIC TRANSLATION ELONGATION	87	0.0000000	0.4428984	0.5154428	0.0000000	0.0000000	0.6795884	0.0512967	0.0000000
1095	TRAFFICKING AND PROCESSING OF ENDOSOMAL TLR	11	0.0008534	0.0789154	0.6536891	0.6504400	0.0001737	0.6584353	0.4064264	0.0064201
534	MET ACTIVATES RAP1 AND RAC1	11	0.0012946	-0.5674576	0.2447823	0.0011181	0.1598323	0.6180021	0.5743404	0.0088796
169	COMPLEX I BIOGENESIS	56	0.0000000	0.2951961	0.5360932	0.0001333	0.0000000	0.6119940	0.1703400	0.0000000
833	RESPONSE OF EIF2AK4 GCN2 TO AMINO ACID DEFICIENCY	94	0.0000000	0.3799739	0.4680230	0.0000000	0.0000000	0.6028480	0.0622601	0.0000000
45	ALPHA LINOLENIC OMEGA3 AND LINOLEIC OMEGA6 ACID METABOLISM	12	0.0022728	0.3720220	0.4719132	0.0256598	0.0046464	0.6009180	0.0706337	0.0143250
1004	SRP DEPENDENT COTRANSLATIONAL PROTEIN TARGETING TO MEMBRANE	106	0.0000000	0.3242635	0.4910307	0.0000000	0.0000000	0.5884369	0.1179223	0.0000000
39	ADP SIGNALLING THROUGH P2Y PURINOCEPTOR 12	19	0.0000308	-0.2223351	0.5447214	0.0934303	0.0000394	0.5883488	0.5423908	0.0004275
443	INITIAL TRIGGERING OF COMPLEMENT	10	0.0085197	0.4741959	0.3371584	0.0094162	0.0648791	0.5818398	0.0969001	0.0395775
123	CD28 DEPENDENT VAV1 PATHWAY	10	0.0045800	-0.4091898	0.4084522	0.0250553	0.0253184	0.5781605	0.5781603	0.0244966
302	EUKARYOTIC TRANSLATION INITIATION	114	0.0000000	0.3056310	0.4747292	0.0000000	0.0000000	0.5646044	0.1195705	0.0000000
655	NUCLEOTIDE LIKE PURINERGIC RECEPTORS	13	0.0024434	0.0355613	0.5555235	0.8243323	0.0005239	0.5566606	0.3676688	0.0149161
830	RESPIRATORY ELECTRON TRANSPORT	102	0.0000000	0.2101871	0.5120165	0.0002465	0.0000000	0.5534794	0.2134256	0.0000000
703	PINK1 PRKN MEDIATED MITOPHAGY	22	0.0000614	0.1877695	0.5206317	0.1274161	0.0000236	0.5534571	0.2353691	0.0007382
1029	SYNTHESIS OF BILE ACIDS AND BILE SALTS VIA 7ALPHA HYDROXYCHOLESTEROL	12	0.0057273	0.3195267	0.4511428	0.0553108	0.0068105	0.5528356	0.0930666	0.0282969
27	ACTIVATION OF THE MRNA UPON BINDING OF THE CAP BINDING COMPLEX AND EIFS AND SUBSEQUENT BINDING TO 43S	59	0.0000000	0.2702378	0.4802650	0.0003316	0.0000000	0.5510744	0.1485117	0.0000000
185	CREB1 PHOSPHORYLATION THROUGH THE ACTIVATION OF ADENYLATE CYCLASE	11	0.0046457	-0.4322950	0.3416882	0.0130441	0.0497435	0.5510262	0.5472888	0.0245107
831	RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CHEMIOSMOTIC COUPLING AND HEAT PRODUCTION BY UNCOUPLING PROTEINS	125	0.0000000	0.2156308	0.5048212	0.0000317	0.0000000	0.5489454	0.2044885	0.0000000
164	COLLAGEN CHAIN TRIMERIZATION	39	0.0000000	0.4231880	-0.3399428	0.0000048	0.0002396	0.5428159	0.5396149	0.0000003
71	ASPARTATE AND ASPARAGINE METABOLISM	10	0.0120759	0.0136397	0.5422249	0.9404699	0.0029861	0.5423964	0.3737662	0.0523441

m2_top <- subset(mres2$enrichment_result,p.adjustMANOVA<0.05)
m2_up <- subset(m2_top,s.dist>0)$set
m2_dn <- subset(m2_top,s.dist<0)$set
#mitch_report(mres2,outfile="mitch2.html",overwrite=TRUE)
write.table(mres2$enrichment_result,file="mitch2.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m2_top_up <- head(subset(m2_top,s.dist>0),10)[,"s.dist"]
names(m2_top_up) <- head(subset(m2_top,s.dist>0),10)[,"set"]
m2_top_dn <- head(subset(m2_top,s.dist<0),10)[,"s.dist"]
names(m2_top_dn) <- head(subset(m2_top,s.dist<0),10)[,"set"]
m2_top_updn <- c(m2_top_up,m2_top_dn)
m2_top_updn <- m2_top_updn[order(m2_top_updn)]

# heatmap
mr2 <- head(mres2$enrichment_result,30)
mr2 <- as.matrix(mr2[,4:5])
rownames(mr2) <- head(mres2$enrichment_result,30)$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mr2,trace="none",col=colfunc(25),scale="none", margins = c(10,28), cexRow=0.7,
  main="Top ranked pathways in amy" , cexCol=0.9 )

Mitch pathway analysis in hippocampus

Mitch for LPS in hip

m3_l <- mitch_import(dge3_l, DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 16927

## Note: no. genes in output = 16911

## Note: estimated proportion of input genes in output = 0.999

mres3_l <- mitch_calc(m3_l, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres3_l$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in hip for LPS") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in hip for LPS
	set	setSize	pANOVA	s.dist	p.adjustANOVA
164	COHESIN LOADING ONTO CHROMATIN	10	0.0035347	-0.5326904	0.0599359
1048	SYNTHESIS SECRETION AND INACTIVATION OF GLUCAGON LIKE PEPTIDE 1 GLP 1	12	0.0014018	-0.5325266	0.0376069
438	INCRETIN SYNTHESIS SECRETION AND INACTIVATION	13	0.0011805	-0.5195426	0.0371346
1057	TERMINATION OF O GLYCAN BIOSYNTHESIS	11	0.0032760	0.5120280	0.0589680
725	POST CHAPERONIN TUBULIN FOLDING PATHWAY	17	0.0002652	0.5108949	0.0173684
191	CROSSLINKING OF COLLAGEN FIBRILS	15	0.0009379	0.4933791	0.0371107
737	PROCESSING AND ACTIVATION OF SUMO	10	0.0082178	-0.4826697	0.1012866
66	APOPTOSIS INDUCED DNA FRAGMENTATION	10	0.0082241	-0.4826223	0.1012866
35	ADENYLATE CYCLASE ACTIVATING PATHWAY	10	0.0106255	0.4665641	0.1106516
201	CYTOSOLIC IRON SULFUR CLUSTER ASSEMBLY	13	0.0037291	0.4645611	0.0623289
1131	TRANSPORT OF CONNEXONS TO THE PLASMA MEMBRANE	13	0.0044157	0.4560121	0.0670955
1047	SYNTHESIS SECRETION AND DEACYLATION OF GHRELIN	12	0.0062768	-0.4556384	0.0863982
125	CD28 DEPENDENT VAV1 PATHWAY	11	0.0114066	-0.4405487	0.1121491
26	ACTIVATION OF SMO	16	0.0028020	0.4315552	0.0555656
579	MITOTIC TELOPHASE CYTOKINESIS	13	0.0112115	-0.4062383	0.1121152
1029	SYNTHESIS OF ACTIVE UBIQUITIN ROLES OF E1 AND E2 ENZYMES	29	0.0001733	-0.4028776	0.0144821
537	MET ACTIVATES RAP1 AND RAC1	11	0.0226278	-0.3969661	0.1647319
826	REPRESSION OF WNT TARGET GENES	14	0.0116740	0.3892830	0.1138216
187	CREB1 PHOSPHORYLATION THROUGH THE ACTIVATION OF ADENYLATE CYCLASE	11	0.0266693	-0.3859387	0.1811960
146	CHK1 CHK2 CDS1 MEDIATED INACTIVATION OF CYCLIN B CDK1 COMPLEX	10	0.0352841	-0.3844506	0.2193881

m3_l_top <- subset(mres3_l$enrichment_result,p.adjustANOVA<0.05)
m3_l_up <- subset(m3_l_top,s.dist>0)$set
m3_l_dn <- subset(m3_l_top,s.dist<0)$set
#mitch_report(mres3_l,outfile="mitch3_l.html",overwrite=TRUE)
write.table(mres3_l$enrichment_result,file="mitch3_l.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m3_l_top_up <- head(subset(m3_l_top,s.dist>0),10)[,"s.dist"]
names(m3_l_top_up) <- head(subset(m3_l_top,s.dist>0),10)[,"set"]
m3_l_top_dn <- head(subset(m3_l_top,s.dist<0),10)[,"s.dist"]
names(m3_l_top_dn) <- head(subset(m3_l_top,s.dist<0),10)[,"set"]
m3_l_top_updn <- c(m3_l_top_up,m3_l_top_dn)
m3_l_top_updn <- m3_l_top_updn[order(m3_l_top_updn)]

par(mar=c(5,25,3,1))
barplot(m3_l_top_updn,horiz=TRUE,las=1,col="darkgray",
  xlab="Enrichment score",cex.names=0.8,xlim=c(-1,1),
  main="Effect of LPS in hip")
grid()

Mitch for OVA in hip

m3_o <- mitch_import(dge2_o, DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 16624

## Note: no. genes in output = 16611

## Note: estimated proportion of input genes in output = 0.999

mres3_o <- mitch_calc(m3_o, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres3_o$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in hip for OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in hip for OVA
	set	setSize	pANOVA	s.dist	p.adjustANOVA
1096	TRAFFICKING AND PROCESSING OF ENDOSOMAL TLR	11	0.0001763	0.6530559	0.0016590
656	NUCLEOTIDE LIKE PURINERGIC RECEPTORS	13	0.0005347	0.5546451	0.0043939
39	ADP SIGNALLING THROUGH P2Y PURINOCEPTOR 12	19	0.0000407	0.5437560	0.0005103
72	ASPARTATE AND ASPARAGINE METABOLISM	10	0.0030231	0.5415336	0.0172096
170	COMPLEX I BIOGENESIS	56	0.0000000	0.5352440	0.0000000
1045	SYNTHESIS OF VERY LONG CHAIN FATTY ACYL COAS	20	0.0000351	0.5344042	0.0004546
1137	TRIGLYCERIDE CATABOLISM	14	0.0006681	0.5251638	0.0053182
704	PINK1 PRKN MEDIATED MITOPHAGY	22	0.0000244	0.5196982	0.0003353
344	FRS MEDIATED FGFR4 SIGNALING	12	0.0019164	0.5173002	0.0121544
921	SHC MEDIATED CASCADE FGFR4	10	0.0047384	0.5157762	0.0243601
302	EUKARYOTIC TRANSLATION ELONGATION	87	0.0000000	0.5144868	0.0000000
188	CRMPS IN SEMA3A SIGNALING	16	0.0003875	-0.5123682	0.0033247
831	RESPIRATORY ELECTRON TRANSPORT	102	0.0000000	0.5111484	0.0000000
153	CITRIC ACID CYCLE TCA CYCLE	22	0.0000354	0.5092805	0.0004546
700	PI 3K CASCADE FGFR4	10	0.0053424	0.5087163	0.0264908
832	RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CHEMIOSMOTIC COUPLING AND HEAT PRODUCTION BY UNCOUPLING PROTEINS	125	0.0000000	0.5039253	0.0000000
718	PLATELET SENSITIZATION BY LDL	15	0.0007660	0.5017755	0.0058008
343	FRS MEDIATED FGFR3 SIGNALING	14	0.0012050	0.4997547	0.0082239
950	SIGNALING BY FGFR3 FUSIONS IN CANCER	10	0.0062828	0.4990422	0.0302976
920	SHC MEDIATED CASCADE FGFR3	12	0.0029553	0.4955419	0.0169061

m3_o_top <- subset(mres3_o$enrichment_result,p.adjustANOVA<0.05)
m3_o_up <- subset(m3_o_top,s.dist>0)$set
m3_o_dn <- subset(m3_o_top,s.dist<0)$set
#mitch_report(mres3_o,outfile="mitch3_o.html",overwrite=TRUE)
write.table(mres3_o$enrichment_result,file="mitch3_o.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m3_o_top_up <- head(subset(m3_o_top,s.dist>0),10)[,"s.dist"]
names(m3_o_top_up) <- head(subset(m3_o_top,s.dist>0),10)[,"set"]
m3_o_top_dn <- head(subset(m3_o_top,s.dist<0),10)[,"s.dist"]
names(m3_o_top_dn) <- head(subset(m3_o_top,s.dist<0),10)[,"set"]
m3_o_top_updn <- c(m3_o_top_up,m3_o_top_dn)
m3_o_top_updn <- m3_o_top_updn[order(m3_o_top_updn)]

par(mar=c(5,25,3,1))
barplot(m3_o_top_updn,horiz=TRUE,las=1,col="darkgray",
  xlab="Enrichment score",cex.names=0.8,xlim=c(-1,1),
  main="Effect of LPS in hip")
grid()

Two dimensional mitch analysis in hip.

ml <- list("LPS"=dge3_l,"OVA"=dge3_o)
m3 <- mitch_import(ml, DEtype="deseq2",geneTable=gt)

## Note: Mean no. genes in input = 16940

## Note: no. genes in output = 16806

## Note: estimated proportion of input genes in output = 0.992

head(m3)

##                      LPS        OVA
## 0610009B22Rik -3.3814425 -2.2697336
## 0610009E02Rik -0.3557659 -0.3519839
## 0610009L18Rik -0.8507498 -0.9139963
## 0610010K14Rik  1.2510787 -0.1751685
## 0610012G03Rik  3.9729326  0.5678859
## 0610030E20Rik  0.1581737  1.4359134

mres3 <- mitch_calc(m3, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(mres3$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in hip for LPS and OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in hip for LPS and OVA
	set	setSize	pMANOVA	s.LPS	s.OVA	p.LPS	p.OVA	s.dist	SD	p.adjustMANOVA
1048	SYNTHESIS SECRETION AND INACTIVATION OF GLUCAGON LIKE PEPTIDE 1 GLP 1	12	0.0004731	-0.5304672	-0.6523262	0.0014631	0.0000911	0.8407883	0.0861673	0.0034594
438	INCRETIN SYNTHESIS SECRETION AND INACTIVATION	13	0.0002930	-0.5173813	-0.6458781	0.0012377	0.0000552	0.8275518	0.0908609	0.0022407
1057	TERMINATION OF O GLYCAN BIOSYNTHESIS	10	0.0133453	0.5267683	0.4883901	0.0039202	0.0074883	0.7183382	0.0271375	0.0536562
1047	SYNTHESIS SECRETION AND DEACYLATION OF GHRELIN	12	0.0065744	-0.4546663	-0.5266464	0.0063889	0.0015832	0.6957571	0.0508976	0.0311421
35	ADENYLATE CYCLASE ACTIVATING PATHWAY	10	0.0172419	0.4668969	0.5136818	0.0105702	0.0049107	0.6941626	0.0330820	0.0652847
66	APOPTOSIS INDUCED DNA FRAGMENTATION	10	0.0224781	-0.4818766	-0.4765063	0.0083240	0.0090747	0.6776897	0.0037974	0.0804263
125	CD28 DEPENDENT VAV1 PATHWAY	11	0.0181447	-0.4388914	-0.4878887	0.0117200	0.0050801	0.6562477	0.0346463	0.0678254
1029	SYNTHESIS OF ACTIVE UBIQUITIN ROLES OF E1 AND E2 ENZYMES	29	0.0000217	-0.4018659	-0.4971729	0.0001799	0.0000036	0.6392785	0.0673922	0.0002283
191	CROSSLINKING OF COLLAGEN FIBRILS	15	0.0041889	0.4934350	0.4061263	0.0009367	0.0064641	0.6390748	0.0617366	0.0214053
371	GENE AND PROTEIN EXPRESSION BY JAK STAT SIGNALING AFTER INTERLEUKIN 12 STIMULATION	31	0.0000000	-0.2950416	-0.5534369	0.0044722	0.0000001	0.6271697	0.1827130	0.0000005
201	CYTOSOLIC IRON SULFUR CLUSTER ASSEMBLY	13	0.0148571	0.4644380	0.3885639	0.0037383	0.0152800	0.6055449	0.0536511	0.0585280
737	PROCESSING AND ACTIVATION OF SUMO	10	0.0295647	-0.4806740	-0.3557037	0.0084871	0.0514563	0.5979738	0.0883673	0.1002628
26	ACTIVATION OF SMO	16	0.0091319	0.4317972	0.4079512	0.0027868	0.0047259	0.5940311	0.0168617	0.0400161
28	ACTIVATION OF THE MRNA UPON BINDING OF THE CAP BINDING COMPLEX AND EIFS AND SUBSEQUENT BINDING TO 43S	59	0.0000000	-0.2071902	-0.5518408	0.0059278	0.0000000	0.5894540	0.2437048	0.0000000
303	EUKARYOTIC TRANSLATION ELONGATION	87	0.0000000	-0.1277581	-0.5729423	0.0395303	0.0000000	0.5870137	0.3147928	0.0000000
189	CRMPS IN SEMA3A SIGNALING	16	0.0021278	0.3208904	0.4886912	0.0262743	0.0007132	0.5846279	0.1186531	0.0125678
910	SEMA3A PLEXIN REPULSION SIGNALING BY INHIBITING INTEGRIN ADHESION	14	0.0150818	0.3652249	0.4470751	0.0179832	0.0037757	0.5772915	0.0578768	0.0590159
304	EUKARYOTIC TRANSLATION INITIATION	114	0.0000000	-0.1675380	-0.5479650	0.0020151	0.0000000	0.5730049	0.2690025	0.0000000
1006	SRP DEPENDENT COTRANSLATIONAL PROTEIN TARGETING TO MEMBRANE	106	0.0000000	-0.1731703	-0.5407016	0.0020803	0.0000000	0.5677554	0.2598839	0.0000000
146	CHK1 CHK2 CDS1 MEDIATED INACTIVATION OF CYCLIN B CDK1 COMPLEX	10	0.0708072	-0.3832103	-0.4122648	0.0358799	0.0239832	0.5628609	0.0205447	0.1922143

m3_top <- subset(mres3$enrichment_result,p.adjustMANOVA<0.05)
m3_up <- subset(m3_top,s.dist>0)$set
m3_dn <- subset(m3_top,s.dist<0)$set
#mitch_report(mres3,outfile="mitch3.html",overwrite=TRUE)
write.table(mres3$enrichment_result,file="mitch3.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m3_top_up <- head(subset(m3_top,s.dist>0),10)[,"s.dist"]
names(m3_top_up) <- head(subset(m3_top,s.dist>0),10)[,"set"]
m3_top_dn <- head(subset(m3_top,s.dist<0),10)[,"s.dist"]
names(m3_top_dn) <- head(subset(m3_top,s.dist<0),10)[,"set"]
m3_top_updn <- c(m3_top_up,m3_top_dn)
m3_top_updn <- m3_top_updn[order(m3_top_updn)]

# heatmap
mr3 <- head(mres3$enrichment_result,30)
mr3 <- as.matrix(mr3[,4:5])
rownames(mr3) <- head(mres3$enrichment_result,30)$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mr3,trace="none",col=colfunc(25),scale="none", margins = c(10,28), cexRow=0.7,
  main="Top ranked pathways in hip" , cexCol=0.9 )

Mitch pathway analysis in PAG

Mitch for LPS in pag

m4_l <- mitch_import(dge4_l, DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 17612

## Note: no. genes in output = 17595

## Note: estimated proportion of input genes in output = 0.999

mres4_l <- mitch_calc(m4_l, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres4_l$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in pag for LPS") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in pag for LPS
	set	setSize	pANOVA	s.dist	p.adjustANOVA
340	FORMATION OF TUBULIN FOLDING INTERMEDIATES BY CCT TRIC	19	0.0000001	0.7073461	0.0000055
47	ALPHA LINOLENIC OMEGA3 AND LINOLEIC OMEGA6 ACID METABOLISM	12	0.0002646	0.6080779	0.0039983
398	GOLGI CISTERNAE PERICENTRIOLAR STACK REORGANIZATION	12	0.0005674	0.5745986	0.0070354
757	PURINE RIBONUCLEOSIDE MONOPHOSPHATE BIOSYNTHESIS	12	0.0007887	0.5595935	0.0091092
860	RHOBTB3 ATPASE CYCLE	10	0.0024461	0.5533125	0.0210319
1138	TRANSPORT OF CONNEXONS TO THE PLASMA MEMBRANE	13	0.0008358	0.5350840	0.0093770
92	BBSOME MEDIATED CARGO TARGETING TO CILIUM	23	0.0000093	0.5339324	0.0002804
154	CITRIC ACID CYCLE TCA CYCLE	22	0.0000249	0.5191590	0.0006362
655	NUCLEOBASE BIOSYNTHESIS	15	0.0006679	0.5073720	0.0080286
505	KERATAN SULFATE BIOSYNTHESIS	24	0.0000254	0.4965141	0.0006362
182	COPI INDEPENDENT GOLGI TO ER RETROGRADE TRAFFIC	44	0.0000000	0.4934192	0.0000012
180	COOPERATION OF PREFOLDIN AND TRIC CCT IN ACTIN AND TUBULIN FOLDING	26	0.0000252	0.4772786	0.0006362
435	HSP90 CHAPERONE CYCLE FOR STEROID HORMONE RECEPTORS SHR	48	0.0000000	0.4651436	0.0000017
507	KERATAN SULFATE KERATIN METABOLISM	30	0.0000146	0.4572199	0.0004093
110	CALNEXIN CALRETICULIN CYCLE	26	0.0000590	0.4550804	0.0011981
391	GLYCOGEN STORAGE DISEASES	12	0.0064361	0.4542551	0.0423560
371	GAP JUNCTION ASSEMBLY	21	0.0003136	0.4542533	0.0046181
284	ENDOSOMAL VACUOLAR PATHWAY	10	0.0131146	-0.4530225	0.0680573
842	RETROGRADE NEUROTROPHIN SIGNALLING	13	0.0047720	0.4520357	0.0339274
392	GLYCOGEN SYNTHESIS	14	0.0036266	0.4490156	0.0273853

m4_l_top <- subset(mres4_l$enrichment_result,p.adjustANOVA<0.05)
m4_l_up <- subset(m4_l_top,s.dist>0)$set
m4_l_dn <- subset(m4_l_top,s.dist<0)$set
#mitch_report(mres4_l,outfile="mitch4_l.html",overwrite=TRUE)
write.table(mres3_l$enrichment_result,file="mitch4_l.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m4_l_top_up <- head(subset(m4_l_top,s.dist>0),10)[,"s.dist"]
names(m4_l_top_up) <- head(subset(m4_l_top,s.dist>0),10)[,"set"]
m4_l_top_dn <- head(subset(m4_l_top,s.dist<0),10)[,"s.dist"]
names(m4_l_top_dn) <- head(subset(m4_l_top,s.dist<0),10)[,"set"]
m4_l_top_updn <- c(m4_l_top_up,m4_l_top_dn)
m4_l_top_updn <- m4_l_top_updn[order(m4_l_top_updn)]

par(mar=c(5,25,3,1))
barplot(m4_l_top_updn,horiz=TRUE,las=1,col="darkgray",
  xlab="Enrichment score",cex.names=0.8,xlim=c(-1,1),
  main="Effect of LPS in pag")
grid()

Mitch for OVA in pag

m4_o <- mitch_import(dge4_o, DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 17662

## Note: no. genes in output = 17643

## Note: estimated proportion of input genes in output = 0.999

mres4_o <- mitch_calc(m4_o, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres4_o$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in pag for OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in pag for OVA
	set	setSize	pANOVA	s.dist	p.adjustANOVA
154	CITRIC ACID CYCLE TCA CYCLE	22	0.0000344	0.5101300	0.0009414
663	OLFACTORY SIGNALING PATHWAY	37	0.0000002	-0.4988533	0.0000111
1046	SYNTHESIS OF PIPS AT THE LATE ENDOSOME MEMBRANE	11	0.0042029	0.4984326	0.0330071
47	ALPHA LINOLENIC OMEGA3 AND LINOLEIC OMEGA6 ACID METABOLISM	12	0.0030441	0.4940162	0.0257983
757	PURINE RIBONUCLEOSIDE MONOPHOSPHATE BIOSYNTHESIS	12	0.0036242	0.4849980	0.0296482
110	CALNEXIN CALRETICULIN CYCLE	26	0.0000235	0.4791045	0.0007681
293	ER QUALITY CONTROL COMPARTMENT ERQC	21	0.0001584	0.4761472	0.0029790
860	RHOBTB3 ATPASE CYCLE	10	0.0093581	0.4745761	0.0581562
104	BRANCHED CHAIN AMINO ACID CATABOLISM	21	0.0001919	0.4700888	0.0033738
1044	SYNTHESIS OF PIPS AT THE EARLY ENDOSOME MEMBRANE	16	0.0011404	0.4697623	0.0126734
448	INHIBITION OF REPLICATION INITIATION OF DAMAGED DNA BY RB1 E2F1	13	0.0039807	0.4612679	0.0316839
74	ASPARTATE AND ASPARAGINE METABOLISM	10	0.0124659	0.4563149	0.0712853
539	MET ACTIVATES RAP1 AND RAC1	10	0.0132563	0.4523224	0.0732017
92	BBSOME MEDIATED CARGO TARGETING TO CILIUM	23	0.0001746	0.4520604	0.0032145
842	RETROGRADE NEUROTROPHIN SIGNALLING	14	0.0048190	0.4351191	0.0363894
111	CAMK IV MEDIATED PHOSPHORYLATION OF CREB	10	0.0187158	0.4293767	0.0896227
907	SCAVENGING OF HEME FROM PLASMA	11	0.0144098	-0.4260745	0.0761198
1051	SYNTHESIS OF VERY LONG CHAIN FATTY ACYL COAS	20	0.0009824	0.4256540	0.0114189
594	N GLYCAN TRIMMING IN THE ER AND CALNEXIN CALRETICULIN CYCLE	35	0.0000156	0.4219186	0.0005813
1073	THE ROLE OF NEF IN HIV 1 REPLICATION AND DISEASE PATHOGENESIS	27	0.0001879	0.4152328	0.0033543

m4_o_top <- subset(mres4_o$enrichment_result,p.adjustANOVA<0.05)
m4_o_up <- subset(m4_o_top,s.dist>0)$set
m4_o_dn <- subset(m4_o_top,s.dist<0)$set
#mitch_report(mres4_o,outfile="mitch4_o.html",overwrite=TRUE)
write.table(mres4_o$enrichment_result,file="mitch4_o.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m4_o_top_up <- head(subset(m4_o_top,s.dist>0),10)[,"s.dist"]
names(m4_o_top_up) <- head(subset(m4_o_top,s.dist>0),10)[,"set"]
m4_o_top_dn <- head(subset(m4_o_top,s.dist<0),10)[,"s.dist"]
names(m4_o_top_dn) <- head(subset(m4_o_top,s.dist<0),10)[,"set"]
m4_o_top_updn <- c(m4_o_top_up,m4_o_top_dn)
m4_o_top_updn <- m4_o_top_updn[order(m4_o_top_updn)]

par(mar=c(5,25,3,1))
barplot(m4_o_top_updn,horiz=TRUE,las=1,col="darkgray",
  xlab="Enrichment score",cex.names=0.8,xlim=c(-1,1),
  main="Effect of LPS in pag")
grid()

Two dimensional mitch analysis in pag.

ml <- list("LPS"=dge4_l,"OVA"=dge4_o)
m4 <- mitch_import(ml, DEtype="deseq2",geneTable=gt)

## Note: Mean no. genes in input = 17637

## Note: no. genes in output = 17476

## Note: estimated proportion of input genes in output = 0.991

head(m4)

##                       LPS         OVA
## 0610009B22Rik  0.02592933  0.45655760
## 0610009E02Rik  0.75195906 -0.40699076
## 0610009L18Rik -1.06444668 -0.18691102
## 0610010K14Rik  0.44760852 -0.37114848
## 0610012G03Rik -0.16331991 -0.34216334
## 0610030E20Rik -0.62496606  0.08425395

mres4 <- mitch_calc(m4, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(mres4$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in pag for LPS and OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in pag for LPS and OVA
	set	setSize	pMANOVA	s.LPS	s.OVA	p.LPS	p.OVA	s.dist	SD	p.adjustMANOVA
47	ALPHA LINOLENIC OMEGA3 AND LINOLEIC OMEGA6 ACID METABOLISM	12	0.0009579	0.6080031	0.4942358	0.0002651	0.0030312	0.7835412	0.0804457	0.0082298
340	FORMATION OF TUBULIN FOLDING INTERMEDIATES BY CCT TRIC	19	0.0000002	0.7072506	0.3203511	0.0000001	0.0156356	0.7764201	0.2735793	0.0000099
756	PURINE RIBONUCLEOSIDE MONOPHOSPHATE BIOSYNTHESIS	12	0.0022717	0.5594175	0.4855417	0.0007918	0.0035867	0.7407419	0.0522381	0.0150215
859	RHOBTB3 ATPASE CYCLE	10	0.0072411	0.5531891	0.4747853	0.0024516	0.0093271	0.7289988	0.0554398	0.0334224
154	CITRIC ACID CYCLE TCA CYCLE	22	0.0000254	0.5190162	0.5105316	0.0000250	0.0000339	0.7280250	0.0059995	0.0004815
92	BBSOME MEDIATED CARGO TARGETING TO CILIUM	23	0.0000282	0.5338437	0.4524026	0.0000093	0.0001727	0.6997551	0.0575875	0.0005266
663	OLFACTORY SIGNALING PATHWAY	37	0.0000001	-0.4715882	-0.4983874	0.0000007	0.0000002	0.6861381	0.0189499	0.0000060
110	CALNEXIN CALRETICULIN CYCLE	26	0.0000326	0.4550584	0.4794578	0.0000590	0.0000232	0.6610279	0.0172530	0.0005901
398	GOLGI CISTERNAE PERICENTRIOLAR STACK REORGANIZATION	12	0.0021591	0.5745629	0.2948255	0.0005678	0.0770110	0.6457899	0.1978042	0.0145958
654	NUCLEOBASE BIOSYNTHESIS	15	0.0027443	0.5072371	0.3844721	0.0006701	0.0099354	0.6364811	0.0868080	0.0171813
841	RETROGRADE NEUROTROPHIN SIGNALLING	13	0.0103429	0.4520371	0.4265678	0.0047719	0.0077447	0.6215284	0.0180095	0.0442675
594	N GLYCAN TRIMMING IN THE ER AND CALNEXIN CALRETICULIN CYCLE	35	0.0000130	0.4196368	0.4223136	0.0000174	0.0000153	0.5953519	0.0018928	0.0002894
104	BRANCHED CHAIN AMINO ACID CATABOLISM	21	0.0008168	0.3588684	0.4703005	0.0044155	0.0001906	0.5915818	0.0787944	0.0073952
293	ER QUALITY CONTROL COMPARTMENT ERQC	21	0.0007289	0.3487853	0.4764660	0.0056603	0.0001568	0.5904837	0.0902839	0.0068629
1050	SYNTHESIS OF VERY LONG CHAIN FATTY ACYL COAS	20	0.0020671	0.3971127	0.4259681	0.0021089	0.0009740	0.5823636	0.0204039	0.0142278
516	LDL CLEARANCE	16	0.0078213	0.4371707	0.3651203	0.0024650	0.0114545	0.5695885	0.0509473	0.0351361
182	COPI INDEPENDENT GOLGI TO ER RETROGRADE TRAFFIC	44	0.0000001	0.4933899	0.2799684	0.0000000	0.0013145	0.5672882	0.1509118	0.0000045
395	GLYCOSPHINGOLIPID METABOLISM	38	0.0000157	0.4065562	0.3944501	0.0000144	0.0000258	0.5664617	0.0085603	0.0003239
776	RECEPTOR MEDIATED MITOPHAGY	11	0.0428258	0.4239180	0.3571663	0.0149128	0.0402577	0.5543232	0.0472006	0.1263308
1048	SYNTHESIS OF PYROPHOSPHATES IN THE CYTOSOL	10	0.0622736	0.4137753	0.3593954	0.0234700	0.0490802	0.5480648	0.0384524	0.1611850

m4_top <- subset(mres4$enrichment_result,p.adjustMANOVA<0.05)
m4_up <- subset(m4_top,s.dist>0)$set
m4_dn <- subset(m4_top,s.dist<0)$set
#mitch_report(mres4,outfile="mitch4.html",overwrite=TRUE)
write.table(mres4$enrichment_result,file="mitch4.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m4_top_up <- head(subset(m4_top,s.dist>0),10)[,"s.dist"]
names(m4_top_up) <- head(subset(m4_top,s.dist>0),10)[,"set"]
m4_top_dn <- head(subset(m4_top,s.dist<0),10)[,"s.dist"]
names(m4_top_dn) <- head(subset(m4_top,s.dist<0),10)[,"set"]
m4_top_updn <- c(m4_top_up,m4_top_dn)
m4_top_updn <- m4_top_updn[order(m4_top_updn)]

# heatmap
mr4 <- head(mres4$enrichment_result,30)
mr4 <- as.matrix(mr4[,4:5])
rownames(mr4) <- head(mres4$enrichment_result,30)$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mr4,trace="none",col=colfunc(25),scale="none", margins = c(10,28), cexRow=0.7,
  main="Top ranked pathways in pag" , cexCol=0.9 )

Mitch pathway analysis in NI

Mitch for LPS in NI

m2_l <- mitch_import(dge2_l, DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 16818

## Note: no. genes in output = 16805

## Note: estimated proportion of input genes in output = 0.999

mres2_l <- mitch_calc(m2_l, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres2_l$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in NI for LPS") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in NI for LPS
	set	setSize	pANOVA	s.dist	p.adjustANOVA
534	MET ACTIVATES RAP1 AND RAC1	11	0.0010392	-0.5710589	0.0295794
1039	SYNTHESIS OF PIPS AT THE LATE ENDOSOME MEMBRANE	11	0.0021453	-0.5344441	0.0510929
24	ACTIVATION OF RAC1	13	0.0017643	-0.5009253	0.0447585
443	INITIAL TRIGGERING OF COMPLEMENT	10	0.0100369	0.4701756	0.1541198
199	CYTOSOLIC IRON SULFUR CLUSTER ASSEMBLY	13	0.0063266	0.4373603	0.1153611
301	EUKARYOTIC TRANSLATION ELONGATION	87	0.0000000	0.4368724	0.0000000
185	CREB1 PHOSPHORYLATION THROUGH THE ACTIVATION OF ADENYLATE CYCLASE	11	0.0121751	-0.4365520	0.1711852
621	NF KB IS ACTIVATED AND SIGNALS SURVIVAL	13	0.0068114	0.4334482	0.1222900
20	ACTIVATION OF IRF3 IRF7 MEDIATED BY TBK1 IKK EPSILON	16	0.0027795	0.4319123	0.0636015
909	SEMA4D INDUCED CELL MIGRATION AND GROWTH CONE COLLAPSE	20	0.0009094	0.4284599	0.0265313
901	SCAVENGING BY CLASS A RECEPTORS	14	0.0075549	0.4123553	0.1306844
917	SEROTONIN RECEPTORS	10	0.0291904	-0.3982971	0.2958961
164	COLLAGEN CHAIN TRIMERIZATION	40	0.0000146	0.3960722	0.0008119
1037	SYNTHESIS OF PIPS AT THE EARLY ENDOSOME MEMBRANE	16	0.0061523	-0.3955938	0.1139641
798	REGULATION OF INNATE IMMUNE RESPONSES TO CYTOSOLIC DNA	13	0.0188749	0.3761223	0.2265723
833	RESPONSE OF EIF2AK4 GCN2 TO AMINO ACID DEFICIENCY	94	0.0000000	0.3738816	0.0000001
45	ALPHA LINOLENIC OMEGA3 AND LINOLEIC OMEGA6 ACID METABOLISM	12	0.0277431	0.3669485	0.2916776
163	COLLAGEN BIOSYNTHESIS AND MODIFYING ENZYMES	63	0.0000010	0.3568954	0.0000704
629	NONSENSE MEDIATED DECAY NMD	109	0.0000000	0.3466105	0.0000001
42	AGGREPHAGY	35	0.0005204	0.3389590	0.0164130

m2_l_top <- subset(mres2_l$enrichment_result,p.adjustANOVA<0.05)
m2_l_up <- subset(m2_l_top,s.dist>0)$set
m2_l_dn <- subset(m2_l_top,s.dist<0)$set
#mitch_report(mres2_l,outfile="mitch2_l.html",overwrite=TRUE)
write.table(mres2_l$enrichment_result,file="mitch2_l.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m2_l_top_up <- head(subset(m2_l_top,s.dist>0),10)[,"s.dist"]
names(m2_l_top_up) <- head(subset(m2_l_top,s.dist>0),10)[,"set"]
m2_l_top_dn <- head(subset(m2_l_top,s.dist<0),10)[,"s.dist"]
names(m2_l_top_dn) <- head(subset(m2_l_top,s.dist<0),10)[,"set"]
m2_l_top_updn <- c(m2_l_top_up,m2_l_top_dn)
m2_l_top_updn <- m2_l_top_updn[order(m2_l_top_updn)]

par(mar=c(5,25,3,1))
barplot(m2_l_top_updn,horiz=TRUE,las=1,col="darkgray",
  xlab="Enrichment score",cex.names=0.8,xlim=c(-1,1),
  main="Effect of LPS in NI")
grid()

Mitch for OVA in NI

m5_o <- mitch_import(dge5_o, DEtype="deseq2",geneTable=gt)

## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.

## Note: Mean no. genes in input = 17391

## Note: no. genes in output = 17372

## Note: estimated proportion of input genes in output = 0.999

mres5_o <- mitch_calc(m5_o, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be
##             statistically significant.

head(mres5_o$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in NI for OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in NI for OVA
	set	setSize	pANOVA	s.dist	p.adjustANOVA
818	REGULATION OF RUNX1 EXPRESSION AND ACTIVITY	17	0.0000127	0.6113207	0.0005565
299	ERYTHROCYTES TAKE UP CARBON DIOXIDE AND RELEASE OXYGEN	10	0.0019843	-0.5647506	0.0220701
638	NOTCH2 INTRACELLULAR DOMAIN REGULATES TRANSCRIPTION	11	0.0039672	0.5016102	0.0377207
1106	TRAFFICKING AND PROCESSING OF ENDOSOMAL TLR	11	0.0046990	-0.4922475	0.0428705
36	ADENYLATE CYCLASE ACTIVATING PATHWAY	10	0.0079378	0.4848059	0.0623909
642	NOTCH4 INTRACELLULAR DOMAIN REGULATES TRANSCRIPTION	20	0.0002383	0.4745908	0.0050174
10	ACETYLCHOLINE REGULATES INSULIN SECRETION	10	0.0095632	0.4732174	0.0711974
1	A TETRASACCHARIDE LINKER SEQUENCE IS REQUIRED FOR GAG SYNTHESIS	26	0.0000690	0.4508865	0.0020855
904	RUNX3 REGULATES NOTCH SIGNALING	13	0.0048864	0.4508191	0.0439775
216	DEFECTIVE B4GALT7 CAUSES EDS PROGEROID TYPE	20	0.0005612	0.4455682	0.0089938
635	NOTCH HLH TRANSCRIPTION PATHWAY	28	0.0000501	0.4427138	0.0016424
792	REGULATION OF COMMISSURAL AXON PATHFINDING BY SLIT AND ROBO	10	0.0170031	0.4358599	0.1022237
640	NOTCH3 INTRACELLULAR DOMAIN REGULATES TRANSCRIPTION	23	0.0003471	0.4308957	0.0064455
571	MITOCHONDRIAL TRANSLATION	93	0.0000000	-0.4303484	0.0000000
1036	SYNTHESIS OF ACTIVE UBIQUITIN ROLES OF E1 AND E2 ENZYMES	29	0.0000619	-0.4297371	0.0019192
679	P75NTR REGULATES AXONOGENESIS	10	0.0196043	-0.4262182	0.1127487
1126	TRANSCRIPTIONAL REGULATION OF PLURIPOTENT STEM CELLS	19	0.0013333	0.4252048	0.0165469
1176	YAP1 AND WWTR1 TAZ STIMULATED GENE EXPRESSION	12	0.0141040	0.4092358	0.0916410
188	CRISTAE FORMATION	31	0.0000821	-0.4086846	0.0023120
798	REGULATION OF GENE EXPRESSION IN LATE STAGE BRANCHING MORPHOGENESIS PANCREATIC BUD PRECURSOR CELLS	15	0.0062783	0.4075551	0.0532530

m5_o_top <- subset(mres5_o$enrichment_result,p.adjustANOVA<0.05)
m5_o_up <- subset(m5_o_top,s.dist>0)$set
m5_o_dn <- subset(m5_o_top,s.dist<0)$set
#mitch_report(mres5_o,outfile="mitch5_o.html",overwrite=TRUE)
write.table(mres5_o$enrichment_result,file="mitch5_o.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m5_o_top_up <- head(subset(m5_o_top,s.dist>0),10)[,"s.dist"]
names(m5_o_top_up) <- head(subset(m5_o_top,s.dist>0),10)[,"set"]
m5_o_top_dn <- head(subset(m5_o_top,s.dist<0),10)[,"s.dist"]
names(m5_o_top_dn) <- head(subset(m5_o_top,s.dist<0),10)[,"set"]
m5_o_top_updn <- c(m5_o_top_up,m5_o_top_dn)
m5_o_top_updn <- m5_o_top_updn[order(m5_o_top_updn)]

par(mar=c(5,25,3,1))
barplot(m5_o_top_updn,horiz=TRUE,las=1,col="darkgray",
  xlab="Enrichment score",cex.names=0.8,xlim=c(-1,1),
  main="Effect of LPS in NI")
grid()

Two dimensional mitch analysis in NI.

ml <- list("LPS"=dge5_l,"OVA"=dge5_o)
m5 <- mitch_import(ml, DEtype="deseq2",geneTable=gt)

## Note: Mean no. genes in input = 17440.5

## Note: no. genes in output = 17243

## Note: estimated proportion of input genes in output = 0.989

head(m5)

##                      LPS        OVA
## 0610009B22Rik -2.3583781 -1.1332185
## 0610009E02Rik -0.4921631 -0.3290659
## 0610009L18Rik -1.5061599  0.1306738
## 0610010K14Rik -0.7608005 -0.3126674
## 0610012G03Rik -0.6314820 -0.2034285
## 0610030E20Rik  1.2199694 -0.5747514

mres5 <- mitch_calc(m5, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(mres5$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in NI for LPS and OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in NI for LPS and OVA
	set	setSize	pMANOVA	s.LPS	s.OVA	p.LPS	p.OVA	s.dist	SD	p.adjustMANOVA
298	ERYTHROCYTES TAKE UP CARBON DIOXIDE AND RELEASE OXYGEN	10	0.0006928	-0.6918354	-0.5629548	0.0001514	0.0020511	0.8919385	0.0911324	0.0056088
35	ADENYLATE CYCLASE ACTIVATING PATHWAY	10	0.0008887	0.6840829	0.4861603	0.0001795	0.0077651	0.8392385	0.1399524	0.0066143
815	REGULATION OF RUNX1 EXPRESSION AND ACTIVITY	17	0.0000342	0.5642497	0.6124736	0.0000562	0.0000123	0.8327674	0.0340995	0.0005431
1103	TRAFFICKING AND PROCESSING OF ENDOSOMAL TLR	11	0.0080603	-0.5150566	-0.4899236	0.0030965	0.0048987	0.7108505	0.0177717	0.0365981
1173	YAP1 AND WWTR1 TAZ STIMULATED GENE EXPRESSION	11	0.0135627	0.4521609	0.4931839	0.0094123	0.0046208	0.6690888	0.0290077	0.0544359
1033	SYNTHESIS OF ACTIVE UBIQUITIN ROLES OF E1 AND E2 ENZYMES	29	0.0000263	-0.4799341	-0.4281359	0.0000077	0.0000659	0.6431463	0.0366269	0.0004614
460	INTERACTION BETWEEN L1 AND ANKYRINS	27	0.0000215	0.5151496	0.3662958	0.0000036	0.0009865	0.6321010	0.1052555	0.0003885
215	DEFECTIVE B4GALT7 CAUSES EDS PROGEROID TYPE	20	0.0013479	0.4258317	0.4471172	0.0009778	0.0005368	0.6174516	0.0150511	0.0093242
677	P75NTR REGULATES AXONOGENESIS	10	0.0397101	-0.4389021	-0.4243254	0.0162480	0.0201542	0.6104811	0.0103073	0.1248638
1	A TETRASACCHARIDE LINKER SEQUENCE IS REQUIRED FOR GAG SYNTHESIS	26	0.0002815	0.3811751	0.4524374	0.0007676	0.0000651	0.5916030	0.0503901	0.0028538
569	MITOCHONDRIAL TRANSLATION	93	0.0000000	-0.4003173	-0.4279645	0.0000000	0.0000000	0.5860098	0.0195496	0.0000000
517	LGI ADAM INTERACTIONS	14	0.0104790	0.4659171	0.3503810	0.0025408	0.0232190	0.5829628	0.0816964	0.0446497
633	NOTCH HLH TRANSCRIPTION PATHWAY	28	0.0002299	0.3609269	0.4439608	0.0009479	0.0000478	0.5721621	0.0587138	0.0024806
187	CRISTAE FORMATION	31	0.0001748	-0.3971760	-0.4064756	0.0001295	0.0000897	0.5683056	0.0065758	0.0019958
965	SIGNALING BY LEPTIN	10	0.0615199	0.4013811	0.4019033	0.0279625	0.0277593	0.5680080	0.0003693	0.1690359
339	FORMATION OF TUBULIN FOLDING INTERMEDIATES BY CCT TRIC	19	0.0059738	-0.3915467	-0.3981348	0.0031303	0.0026613	0.5584086	0.0046585	0.0292814
158	CLASS I PEROXISOMAL MEMBRANE PROTEIN IMPORT	20	0.0027202	-0.4438774	-0.3317250	0.0005891	0.0102258	0.5541377	0.0793037	0.0164895
901	RUNX3 REGULATES NOTCH SIGNALING	13	0.0183734	0.3160230	0.4522613	0.0485198	0.0047513	0.5517344	0.0963350	0.0703814
640	NOTCH4 INTRACELLULAR DOMAIN REGULATES TRANSCRIPTION	20	0.0006859	0.2624862	0.4759740	0.0421532	0.0002285	0.5435534	0.1509587	0.0056018
636	NOTCH2 INTRACELLULAR DOMAIN REGULATES TRANSCRIPTION	11	0.0058700	0.2054634	0.5030704	0.2380562	0.0038630	0.5434106	0.2104399	0.0290049

m5_top <- subset(mres5$enrichment_result,p.adjustMANOVA<0.05)
m5_up <- subset(m5_top,s.dist>0)$set
m5_dn <- subset(m5_top,s.dist<0)$set
#mitch_report(mres5,outfile="mitch5.html",overwrite=TRUE)
write.table(mres5$enrichment_result,file="mitch5.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m5_top_up <- head(subset(m5_top,s.dist>0),10)[,"s.dist"]
names(m5_top_up) <- head(subset(m5_top,s.dist>0),10)[,"set"]
m5_top_dn <- head(subset(m5_top,s.dist<0),10)[,"s.dist"]
names(m5_top_dn) <- head(subset(m5_top,s.dist<0),10)[,"set"]
m5_top_updn <- c(m5_top_up,m5_top_dn)
m5_top_updn <- m5_top_updn[order(m5_top_updn)]

# heatmap
mr5 <- head(mres5$enrichment_result,30)
mr5 <- as.matrix(mr5[,4:5])
rownames(mr5) <- head(mres5$enrichment_result,30)$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mr5,trace="none",col=colfunc(25),scale="none", margins = c(10,28), cexRow=0.7,
  main="Top ranked pathways in NI" , cexCol=0.9 )

Mitch pathway analysis - LPS all tissues

ml <- list("hyp"=dge1_l,
  "amy"=dge2_l,
  "hip"=dge3_l, 
  "pag"=dge4_l, 
  "ni"=dge5_l)

m6 <- mitch_import(ml, DEtype="deseq2",geneTable=gt)

## Note: Mean no. genes in input = 17188.4

## Note: no. genes in output = 15915

## Note: estimated proportion of input genes in output = 0.926

head(m6)

##                      hyp        amy        hip         pag         ni
## 0610009B22Rik -0.3304181  0.2684894 -3.3814425  0.02592933 -2.3583781
## 0610009E02Rik  0.7596333 -0.7521235 -0.3557659  0.75195906 -0.4921631
## 0610009L18Rik -1.2521412 -1.1838923 -0.8507498 -1.06444668 -1.5061599
## 0610010K14Rik  0.4193627 -0.0433030  1.2510787  0.44760852 -0.7608005
## 0610012G03Rik  1.4462155 -0.1930476  3.9729326 -0.16331991 -0.6314820
## 0610030E20Rik -0.2426628  1.0496464  0.1581737 -0.62496606  1.2199694

mres6 <- mitch_calc(m6, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(mres6$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in all tissues caused by LPS") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in all tissues caused by LPS
	set	setSize	pMANOVA	s.hyp	s.amy	s.hip	s.pag	s.ni	p.hyp	p.amy	p.hip	p.pag	p.ni	s.dist	SD	p.adjustMANOVA
332	FORMATION OF TUBULIN FOLDING INTERMEDIATES BY CCT TRIC	19	0.0000000	0.5822782	0.1560472	0.1543718	0.6954017	-0.3830159	0.0000111	0.2390518	0.2441327	0.0000002	0.0038508	1.0087190	0.4263417	0.0000000
1118	TRANSPORT OF CONNEXONS TO THE PLASMA MEMBRANE	13	0.0000173	0.5932684	0.3166027	0.4487196	0.5229724	-0.1371864	0.0002123	0.0481161	0.0050908	0.0010950	0.3918199	0.9725606	0.2903895	0.0001283
713	POST CHAPERONIN TUBULIN FOLDING PATHWAY	17	0.0000011	0.6303420	0.3164734	0.5029341	0.3588835	-0.0844797	0.0000068	0.0238916	0.0003304	0.0104181	0.5465488	0.9414687	0.2701030	0.0000104
33	ADENYLATE CYCLASE ACTIVATING PATHWAY	10	0.0001919	-0.0883999	0.2081232	0.4566237	0.1617479	0.6784785	0.6283875	0.2544891	0.0124090	0.3758371	0.0002028	0.8637884	0.2935960	0.0011038
187	CROSSLINKING OF COLLAGEN FIBRILS	13	0.0177579	0.3350715	0.4012945	0.5670985	0.0343934	0.3111268	0.0364702	0.0122413	0.0003994	0.8300153	0.0521183	0.8324023	0.1930438	0.0515781
529	MET ACTIVATES RAP1 AND RAC1	10	0.0644943	-0.3635586	-0.5403458	-0.3609305	0.1324992	-0.3217227	0.0465211	0.0030880	0.0481273	0.4681766	0.0781467	0.8218764	0.2513084	0.1390960
894	SEALING OF THE NUCLEAR ENVELOPE NE BY ESCRT III	23	0.0000005	0.5135808	0.2584401	0.3150505	0.4227394	-0.1673087	0.0000201	0.0319369	0.0089178	0.0004491	0.1649244	0.7978183	0.2627031	0.0000052
121	CD28 DEPENDENT VAV1 PATHWAY	11	0.0259587	-0.1266234	-0.3169586	-0.4383688	0.3942290	-0.3995679	0.4671687	0.0687421	0.0118218	0.0235806	0.0217584	0.7897686	0.3414227	0.0689846
197	CYTOSOLIC IRON SULFUR CLUSTER ASSEMBLY	13	0.0137116	0.3881757	0.4365392	0.4570591	-0.1598735	0.1895649	0.0153845	0.0064267	0.0043267	0.3183090	0.2366974	0.7820769	0.2586774	0.0419327
23	ACTIVATION OF RAC1	13	0.0213128	-0.2225555	-0.5009723	-0.3742248	0.2292019	-0.2882173	0.1647648	0.0017628	0.0194863	0.1525120	0.0719979	0.7590462	0.2777268	0.0596552
183	CREB1 PHOSPHORYLATION THROUGH THE ACTIVATION OF ADENYLATE CYCLASE	11	0.0165125	-0.3768635	-0.4367798	-0.3853119	0.2699378	-0.1380785	0.0304523	0.0121316	0.0269197	0.1211234	0.4278525	0.7570997	0.2938798	0.0488194
362	GAP JUNCTION ASSEMBLY	21	0.0000217	0.4383385	0.2261051	0.3194856	0.4386022	-0.1332099	0.0005065	0.0729060	0.0112697	0.0005026	0.2907131	0.7452865	0.2361048	0.0001555
1017	SYNTHESIS OF ACTIVE UBIQUITIN ROLES OF E1 AND E2 ENZYMES	29	0.0000107	-0.0486831	-0.2218088	-0.4020673	0.3228694	-0.4739632	0.6500898	0.0387355	0.0001787	0.0026213	0.0000100	0.7362833	0.3187602	0.0000844
176	COOPERATION OF PREFOLDIN AND TRIC CCT IN ACTIN AND TUBULIN FOLDING	26	0.0000000	0.3865277	0.1647439	0.0234076	0.4606525	-0.3841797	0.0006466	0.1460230	0.8363651	0.0000479	0.0006975	0.7327256	0.3361871	0.0000003
44	ALPHA LINOLENIC OMEGA3 AND LINOLEIC OMEGA6 ACID METABOLISM	12	0.0022542	0.1522040	0.3666185	0.0550630	0.5903813	0.0983462	0.3613373	0.0278868	0.7412277	0.0003981	0.5553169	0.7202979	0.2236017	0.0092404
139	CGMP EFFECTS	15	0.0037150	-0.5594549	-0.3014927	-0.2268679	-0.2409392	0.0179958	0.0001756	0.0432289	0.1282361	0.1062140	0.9039657	0.7167512	0.2062247	0.0142203
772	REDUCTION OF CYTOSOLIC CA LEVELS	11	0.0185558	-0.5540287	-0.3044288	-0.1996296	-0.2657079	0.0577099	0.0014636	0.0804376	0.2516620	0.1270637	0.7403589	0.7165248	0.2195527	0.0531030
1036	SYNTHESIS SECRETION AND INACTIVATION OF GLUCAGON LIKE PEPTIDE 1 GLP 1	12	0.0194825	-0.1096124	-0.1519105	-0.5265254	0.0254355	-0.4450104	0.5109369	0.3622638	0.0015874	0.8787585	0.0076045	0.7148440	0.2341640	0.0553359
680	PHASE 2 PLATEAU PHASE	10	0.1069460	0.1580509	0.1611443	0.3951713	-0.1846589	0.5064948	0.3868463	0.3776212	0.0304834	0.3119985	0.0055468	0.7055097	0.2659891	0.1984424
725	PROCESSING AND ACTIVATION OF SUMO	10	0.0315656	0.0187111	-0.2962087	-0.4767180	0.2294499	-0.3573342	0.9184037	0.1048393	0.0090451	0.2090122	0.0504005	0.7040488	0.2915753	0.0796721

m6_top <- subset(mres6$enrichment_result,p.adjustMANOVA<0.05)
m6_up <- subset(m6_top,s.dist>0)$set
m6_dn <- subset(m6_top,s.dist<0)$set
#mitch_report(mres6,outfile="mitch6.html",overwrite=TRUE)
write.table(mres6$enrichment_result,file="mitch6.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m6_top_up <- head(subset(m6_top,s.dist>0),10)[,"s.dist"]
names(m6_top_up) <- head(subset(m6_top,s.dist>0),10)[,"set"]
m6_top_dn <- head(subset(m6_top,s.dist<0),10)[,"s.dist"]
names(m6_top_dn) <- head(subset(m6_top,s.dist<0),10)[,"set"]
m6_top_updn <- c(m6_top_up,m6_top_dn)
m6_top_updn <- m6_top_updn[order(m6_top_updn)]

# heatmap
mr6 <- head(mres6$enrichment_result,30)
mr6 <- as.matrix(mr6[,4:8])
rownames(mr6) <- head(mres6$enrichment_result,30)$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mr6,trace="none",col=colfunc(25),scale="none", margins = c(10,28), cexRow=0.7,
  main="Top ranked pathways LPS all tissues" , cexCol=0.9 )

Mitch pathway analysis - OVA all tissues

ml <- list("hyp"=dge1_o,
  "amy"=dge2_o,
  "hip"=dge3_o,
  "pag"=dge4_o,
  "ni"=dge5_o)

m7 <- mitch_import(ml, DEtype="deseq2",geneTable=gt)

## Note: Mean no. genes in input = 17148.8

## Note: no. genes in output = 15857

## Note: estimated proportion of input genes in output = 0.925

head(m7)

##                       hyp          amy        hip         pag         ni
## 0610009B22Rik  0.09722066  1.831236577 -2.2697336  0.45655760 -1.1332185
## 0610009E02Rik  1.26326253 -0.103702849 -0.3519839 -0.40699076 -0.3290659
## 0610009L18Rik -0.66374969 -0.075637695 -0.9139963 -0.18691102  0.1306738
## 0610010K14Rik -0.45368981  0.417670579 -0.1751685 -0.37114848 -0.3126674
## 0610012G03Rik  0.12236243  0.008709563  0.5678859 -0.34216334 -0.2034285
## 0610030E20Rik -0.04562294 -1.609523779  1.4359134  0.08425395 -0.5747514

mres7 <- mitch_calc(m7, genesets, priority="effect")

## Note: Enrichments with large effect sizes may not be 
##             statistically significant.

head(mres7$enrichment_result,20) %>%
  kbl(caption = "Top gene pathway differences in all tissues caused by OVA") %>%
  kable_paper("hover", full_width = F)

Top gene pathway differences in all tissues caused by OVA
	set	setSize	pMANOVA	s.hyp	s.amy	s.hip	s.pag	s.ni	p.hyp	p.amy	p.hip	p.pag	p.ni	s.dist	SD	p.adjustMANOVA
301	EUKARYOTIC TRANSLATION ELONGATION	87	0.0000000	0.7588262	0.5064512	-0.5678613	-0.3418057	-0.0059097	0.0000000	0.0000000	0.0000000	0.0000000	0.9241530	1.1276713	0.5583871	0.0000000
998	SRP DEPENDENT COTRANSLATIONAL PROTEIN TARGETING TO MEMBRANE	106	0.0000000	0.6787601	0.4821449	-0.5361121	-0.2519301	-0.0297124	0.0000000	0.0000000	0.0000000	0.0000075	0.5974752	1.0222264	0.5053208	0.0000000
828	RESPONSE OF EIF2AK4 GCN2 TO AMINO ACID DEFICIENCY	94	0.0000000	0.6580641	0.4592548	-0.5303721	-0.3405349	-0.0378344	0.0000000	0.0000000	0.0000000	0.0000000	0.5264627	1.0211040	0.5084172	0.0000000
169	COMPLEX I BIOGENESIS	56	0.0000000	0.5240378	0.5286035	-0.4897972	-0.3222988	-0.3638965	0.0000000	0.0000000	0.0000000	0.0000303	0.0000025	1.0150053	0.5067526	0.0000000
302	EUKARYOTIC TRANSLATION INITIATION	114	0.0000000	0.6323601	0.4657442	-0.5438585	-0.2323617	-0.0824059	0.0000000	0.0000000	0.0000000	0.0000185	0.1289137	0.9865910	0.4903804	0.0000000
325	FORMATION OF ATP BY CHEMIOSMOTIC COUPLING	18	0.0000007	0.5445420	0.4758367	-0.4187764	-0.2660101	-0.3879243	0.0000634	0.0004739	0.0020991	0.0507461	0.0043833	0.9589405	0.4793274	0.0000086
27	ACTIVATION OF THE MRNA UPON BINDING OF THE CAP BINDING COMPLEX AND EIFS AND SUBSEQUENT BINDING TO 43S	59	0.0000000	0.5883366	0.4716473	-0.5476128	-0.1011306	-0.1302310	0.0000000	0.0000000	0.0000000	0.1793102	0.0837500	0.9463924	0.4690057	0.0000000
825	RESPIRATORY ELECTRON TRANSPORT	102	0.0000000	0.4632093	0.5044101	-0.4530750	-0.2710089	-0.3775272	0.0000000	0.0000000	0.0000000	0.0000023	0.0000000	0.9435263	0.4708108	0.0000000
782	REGULATION OF COMMISSURAL AXON PATHFINDING BY SLIT AND ROBO	10	0.0011531	-0.6278665	-0.3110115	0.4449927	0.0985549	0.4317032	0.0005855	0.0885893	0.0148264	0.5894748	0.0180880	0.9407676	0.4703135	0.0056676
1089	TRAFFICKING AND PROCESSING OF ENDOSOMAL TLR	11	0.0032723	0.0864113	0.6469886	-0.4574369	0.0249791	-0.4885202	0.6197636	0.0002025	0.0086156	0.8859484	0.0050240	0.9351928	0.4657114	0.0131346
826	RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CHEMIOSMOTIC COUPLING AND HEAT PRODUCTION BY UNCOUPLING PROTEINS	125	0.0000000	0.4662110	0.4971360	-0.4333618	-0.2414320	-0.3586402	0.0000000	0.0000000	0.0000000	0.0000032	0.0000000	0.9160849	0.4577743	0.0000000
72	ASPARTATE AND ASPARAGINE METABOLISM	10	0.0010542	-0.3786837	0.5351802	-0.1915441	0.4258093	-0.3953808	0.0381296	0.0033837	0.2943000	0.0197252	0.0303956	0.8967429	0.4483703	0.0052485
900	SELENOAMINO ACID METABOLISM	107	0.0000000	0.5893108	0.4098751	-0.4925596	-0.2055054	-0.0288895	0.0000000	0.0000000	0.0000000	0.0002426	0.6059872	0.8949675	0.4433240	0.0000000
806	REGULATION OF RUNX1 EXPRESSION AND ACTIVITY	17	0.0004877	0.0136141	-0.3894014	0.4431670	-0.2126634	0.6067068	0.9225960	0.0054441	0.0015594	0.1290543	0.0000148	0.8726595	0.4239553	0.0026561
1021	SYNTHESIS OF ACTIVE UBIQUITIN ROLES OF E1 AND E2 ENZYMES	29	0.0000045	0.2194060	0.4861747	-0.4948977	0.1091475	-0.4279670	0.0408958	0.0000059	0.0000040	0.3091226	0.0000664	0.8511742	0.4248996	0.0000447
152	CITRIC ACID CYCLE TCA CYCLE	22	0.0000256	0.0600683	0.5015070	-0.3268996	0.4826076	-0.3581652	0.6258176	0.0000466	0.0079547	0.0000890	0.0036388	0.8503959	0.4175464	0.0002169
1158	YAP1 AND WWTR1 TAZ STIMULATED GENE EXPRESSION	10	0.0027124	0.6407648	0.0984540	0.2451316	0.0291285	0.4779327	0.0004499	0.5898562	0.1795474	0.8732901	0.0088707	0.8423956	0.2572784	0.0113589
186	CRISTAE FORMATION	31	0.0000001	0.3863752	0.4229585	-0.4352128	-0.1599940	-0.4056004	0.0001969	0.0000459	0.0000274	0.1232308	0.0000930	0.8412489	0.4184397	0.0000007
625	NONSENSE MEDIATED DECAY NMD	109	0.0000000	0.5291588	0.3640689	-0.4406641	-0.2886506	-0.0092885	0.0000000	0.0000000	0.0000000	0.0000002	0.8670760	0.8307500	0.4139335	0.0000000
1109	TRANSCRIPTIONAL REGULATION OF PLURIPOTENT STEM CELLS	15	0.0000158	0.3654757	0.4781214	0.0412238	0.3560662	0.4261246	0.0142632	0.0013455	0.7822523	0.0169652	0.0042723	0.8199008	0.1706251	0.0001422

m7_top <- subset(mres7$enrichment_result,p.adjustMANOVA<0.05)
m7_up <- subset(m7_top,s.dist>0)$set
m7_dn <- subset(m7_top,s.dist<0)$set
#mitch_report(mres7,outfile="mitch7.html",overwrite=TRUE)
write.table(mres7$enrichment_result,file="mitch7.tsv",quote=FALSE,sep="\t",row.names=FALSE)

m7_top_up <- head(subset(m7_top,s.dist>0),10)[,"s.dist"]
names(m7_top_up) <- head(subset(m7_top,s.dist>0),10)[,"set"]
m7_top_dn <- head(subset(m7_top,s.dist<0),10)[,"s.dist"]
names(m7_top_dn) <- head(subset(m7_top,s.dist<0),10)[,"set"]
m7_top_updn <- c(m7_top_up,m7_top_dn)
m7_top_updn <- m7_top_updn[order(m7_top_updn)]

# heatmap
mr7 <- head(mres7$enrichment_result,30)
mr7 <- as.matrix(mr7[,4:8])
rownames(mr7) <- head(mres7$enrichment_result,30)$set

colfunc <- colorRampPalette(c("blue", "white", "red"))

heatmap.2(mr7,trace="none",col=colfunc(25),scale="none", margins = c(10,28), cexRow=0.7,
  main="Top ranked pathways OVA all tissues" , cexCol=0.9 )

Conclusion

There was some variability among samples in the hypothalamus, in particular some OVA samples could be classified as outliers (EA9-1 and EA3-1). Further investigation would be required to understand the origin of this variability. LPS resulted in the differential expression of just 30 genes in the hypothalamus. OVA was variable, so no significant genes were identified. Analysis with mitch identified some interesting pathways as seen in the mitch heatmap (see the mitch reports for details). The mitch software can perform multi-contrast pathway analysis and will be important to contrast the effect of the treatments in different tissues.

References

Bibliography

 Babraham bioinformatics - FastQC A quality control tool for high throughput sequence data. Babraham.ac.uk, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed February 23, 2022).

 Frankish A, Diekhans M, Jungreis I, et al. GENCODE 2021. Nucleic Acids Res 2021; 49: D916–D923.

 Jiang H, Lei R, Ding S-W, et al. Skewer: a fast and accurate adapter trimmer for next-generation sequencing paired-end reads. BMC Bioinformatics 2014; 15: 182.

 Bray NL, Pimentel H, Melsted P, et al. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol 2016; 34: 525–527.

 Ewels P, Magnusson M, Lundin S, et al. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 2016; 32: 3047–3048.

 Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014; 15: 550.

 Liberzon A, Birger C, Thorvaldsdóttir H, et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 2015; 1: 417–425.

 Jassal B, Matthews L, Viteri G, et al. The reactome pathway knowledgebase. Nucleic Acids Res 2020; 48: D498–D503.

 Dolgalev I. MSigDB Gene Sets for Multiple Organisms in a Tidy Data Format [R package msigdbr version 7.4.1], https://cran.r-project.org/web/packages/msigdbr/index.html (2021, accessed February 23, 2022).

Kaspi A, Ziemann M. Mitch: Multi-contrast pathway enrichment for multi-omics and single-cell profiling data. BMC Genomics 2020; 21: 447.

Session information

For reproducibility.

sessionInfo()

## R version 4.1.3 (2022-03-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.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       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] beeswarm_0.4.0              vioplot_0.3.7              
##  [3] sm_2.2-5.7                  kableExtra_1.3.4           
##  [5] topconfects_1.8.0           limma_3.48.3               
##  [7] eulerr_6.1.1                mitch_1.4.1                
##  [9] MASS_7.3-55                 fgsea_1.18.0               
## [11] gplots_3.1.1                DESeq2_1.32.0              
## [13] SummarizedExperiment_1.22.0 Biobase_2.52.0             
## [15] MatrixGenerics_1.4.3        matrixStats_0.61.0         
## [17] GenomicRanges_1.44.0        GenomeInfoDb_1.28.4        
## [19] IRanges_2.26.0              S4Vectors_0.30.2           
## [21] BiocGenerics_0.38.0         reshape2_1.4.4             
## [23] forcats_0.5.1               stringr_1.4.0              
## [25] dplyr_1.0.8                 purrr_0.3.4                
## [27] readr_2.1.2                 tidyr_1.2.0                
## [29] tibble_3.1.6                ggplot2_3.3.5              
## [31] tidyverse_1.3.1             zoo_1.8-9                  
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1           backports_1.4.1        fastmatch_1.1-3       
##   [4] systemfonts_1.0.3      plyr_1.8.6             polylabelr_0.2.0      
##   [7] splines_4.1.3          BiocParallel_1.26.2    digest_0.6.29         
##  [10] htmltools_0.5.2        fansi_1.0.2            magrittr_2.0.2        
##  [13] memoise_2.0.1          tzdb_0.2.0             Biostrings_2.60.2     
##  [16] annotate_1.70.0        modelr_0.1.8           svglite_2.1.0         
##  [19] colorspace_2.0-2       blob_1.2.2             rvest_1.0.2           
##  [22] haven_2.4.3            xfun_0.29              crayon_1.4.2          
##  [25] RCurl_1.98-1.6         jsonlite_1.7.3         genefilter_1.74.1     
##  [28] survival_3.2-13        glue_1.6.1             polyclip_1.10-0       
##  [31] gtable_0.3.0           zlibbioc_1.38.0        XVector_0.32.0        
##  [34] webshot_0.5.2          DelayedArray_0.18.0    scales_1.1.1          
##  [37] DBI_1.1.2              GGally_2.1.2           Rcpp_1.0.8            
##  [40] viridisLite_0.4.0      xtable_1.8-4           bit_4.0.4             
##  [43] htmlwidgets_1.5.4      httr_1.4.2             RColorBrewer_1.1-2    
##  [46] ellipsis_0.3.2         pkgconfig_2.0.3        reshape_0.8.8         
##  [49] XML_3.99-0.8           sass_0.4.0             dbplyr_2.1.1          
##  [52] locfit_1.5-9.4         utf8_1.2.2             tidyselect_1.1.1      
##  [55] rlang_1.0.1            later_1.3.0            AnnotationDbi_1.54.1  
##  [58] munsell_0.5.0          cellranger_1.1.0       tools_4.1.3           
##  [61] cachem_1.0.6           cli_3.1.1              generics_0.1.2        
##  [64] RSQLite_2.2.9          broom_0.7.12           evaluate_0.14         
##  [67] fastmap_1.1.0          yaml_2.2.2             knitr_1.37            
##  [70] bit64_4.0.5            fs_1.5.2               caTools_1.18.2        
##  [73] KEGGREST_1.32.0        mime_0.12              xml2_1.3.3            
##  [76] compiler_4.1.3         rstudioapi_0.13        png_0.1-7             
##  [79] reprex_2.0.1           geneplotter_1.70.0     bslib_0.3.1           
##  [82] stringi_1.7.6          highr_0.9              lattice_0.20-45       
##  [85] Matrix_1.4-0           vctrs_0.3.8            pillar_1.7.0          
##  [88] lifecycle_1.0.1        jquerylib_0.1.4        data.table_1.14.2     
##  [91] bitops_1.0-7           httpuv_1.6.5           R6_2.5.1              
##  [94] promises_1.2.0.1       KernSmooth_2.23-20     echarts4r_0.4.3       
##  [97] gridExtra_2.3          gtools_3.9.2           assertthat_0.2.1      
## [100] withr_2.4.3            GenomeInfoDbData_1.2.6 hms_1.1.1             
## [103] grid_4.1.3             rmarkdown_2.11         shiny_1.7.1           
## [106] lubridate_1.8.0