Methods

Rats were kept in sedentary conditions or were trained. RNA was isolated from whole cell and mito fractions. Libraries were generated by the Deakin Genomics Facility and sequenced on the NovaSeq system.

Fastqc and MultiQC were run to summarise the QC checks that were done.

Reads were then were mapped to the rat genome build mRatBN7.2 with Ensembl v106 annotations using STAR aligner version 2.7.1a.

imported to R for analysis with DESeq2. Pathway level analysis was then done using mitch with Reactome gene sets.

Read counts

Import the STAR transcript counts. We can also include some info out of the Ensembl GTF file including gene name and gene class.

# libraries
library("reshape2")
library("DESeq2")
library("mitch")
library("gplots")
library("kableExtra")

# import the 3 column table
tmp<-read.table("3col.tsv.gz",header=F)

# convert the 3 col table into a standard count matrix 
x<-as.matrix(acast(tmp, V2~V1, value.var="V3"))
# tidy up the column headers
colnames(x)<-sapply(strsplit(colnames(x),"_"),"[[",1)
head(x)
##                      1  10  11  12  13  14  15  16  17  18  19   2  20  21  22
## ENSRNOG00000000001  18  32  15  16  31  18  12  15  49  20  20  22  16  16  18
## ENSRNOG00000000007   8   5  10   6  12   1   6   3   9   8  19   3   5  21   5
## ENSRNOG00000000008   6  46  46  16  25  44  39  23  73  40  54  21  30  64  11
## ENSRNOG00000000009   1   1   0   0   0   1   0   0   0   0   0   0   0   2   0
## ENSRNOG00000000010   0   0   1   1   0   1   0   1   1   0   1   4   2   1   6
## ENSRNOG00000000012 278 678 336 527 702 446 307 347 368 482 501 322 319 906 598
##                     23  24  25  26  27  28  29   3  30  31  32 33 34 35 36 37
## ENSRNOG00000000001  24  29  20  12  27  16  16  16  10  24  26 14 12 30 13 28
## ENSRNOG00000000007   8   2  12   5  11  10   5  11   9   9   6 34 21 41 34 37
## ENSRNOG00000000008  38  14  35  44  62  29  31  32  34  28  22 14 10 17  6 17
## ENSRNOG00000000009   0   0   0   0   0   0   0   0   0   0   0 18  7 20 11 23
## ENSRNOG00000000010   3   3   0   0   0   2   3   1   0   4   0 23 11 24 16 30
## ENSRNOG00000000012 834 630 619 648 867 758 472 493 296 468 289 19 11 10  6  5
##                    38 39   4 40 41 42 43 44 45 46 47 48 49   5 50 51 52 53 54
## ENSRNOG00000000001 35 22  13 24 16 17 29 27 36 14  8 26 25  19 19 37 31 31 20
## ENSRNOG00000000007 42 28   6 17 25 18 44 35 41 38 41 32 40   8 49 52 48 33 42
## ENSRNOG00000000008 25 13  18 16  6  8 40 15 17  9 31  6 19  25 17 25 18 13 14
## ENSRNOG00000000009 29 26   0 22  5  8 20 10 22 13 21  7 22   0 27 36 33 22 33
## ENSRNOG00000000010 18 12   0 21 17 14 24 13 20 25 21 13 18   0 25 37 16 24 23
## ENSRNOG00000000012 18  2 247  7  4  6 19 14 13 10  6 11  4 832 15  8 15 10 15
##                    55 56 57 58 59   6 60 61 62 63 64   7   8   9 95 96 test
## ENSRNOG00000000001 14 16 21 26 14  18 24 32 33 30 30  14  15  13  0  0    0
## ENSRNOG00000000007 26 38 58 38 34  19 49 44 53 48 41   4   3   6  0  0    0
## ENSRNOG00000000008 25  5 20 12 12  28 22 27 15 18 21  14  14  25  0  0    0
## ENSRNOG00000000009 24 17 28 15 15   0 10 23 21 25 18   0   0   0  0  0    0
## ENSRNOG00000000010 24 24 17 27 27   0 25 27 26 30 32   0   2   0  0  0    0
## ENSRNOG00000000012  6  7 12  8  5 646 10 10  9 10  4 351 269 438  3 33    0
#dont forget gene names
g<-read.table("../ref2/Rattus_norvegicus.mRatBN7.2.106.gnames.txt",row.names=1)
g$gene_ID <- paste(rownames(g),g$V2,g$V3)
head(g)
##                          V2             V3
## ENSRNOG00000066169       na protein_coding
## ENSRNOG00000070168    Olr56 protein_coding
## ENSRNOG00000070901     Irgq protein_coding
## ENSRNOG00000018029    Doc2g protein_coding
## ENSRNOG00000031391 Ceacam16 protein_coding
## ENSRNOG00000055342       U7          snRNA
##                                                       gene_ID
## ENSRNOG00000066169       ENSRNOG00000066169 na protein_coding
## ENSRNOG00000070168    ENSRNOG00000070168 Olr56 protein_coding
## ENSRNOG00000070901     ENSRNOG00000070901 Irgq protein_coding
## ENSRNOG00000018029    ENSRNOG00000018029 Doc2g protein_coding
## ENSRNOG00000031391 ENSRNOG00000031391 Ceacam16 protein_coding
## ENSRNOG00000055342                ENSRNOG00000055342 U7 snRNA
g <- g[,ncol(g),drop=FALSE] 

x<-merge(g,x,by=0)

x[,1]=NULL
rownames(x) <- x[,1]
x[,1]=NULL
x$test=NULL
head(x)
##                                            1  10  11  12  13  14  15  16  17
## ENSRNOG00000000001 Arsj protein_coding    18  32  15  16  31  18  12  15  49
## ENSRNOG00000000007 Gad1 protein_coding     8   5  10   6  12   1   6   3   9
## ENSRNOG00000000008 Alx4 protein_coding     6  46  46  16  25  44  39  23  73
## ENSRNOG00000000009 Tmco5b protein_coding   1   1   0   0   0   1   0   0   0
## ENSRNOG00000000010 Cbln1 protein_coding    0   0   1   1   0   1   0   1   1
## ENSRNOG00000000012 Tcf15 protein_coding  278 678 336 527 702 446 307 347 368
##                                           18  19   2  20  21  22  23  24  25
## ENSRNOG00000000001 Arsj protein_coding    20  20  22  16  16  18  24  29  20
## ENSRNOG00000000007 Gad1 protein_coding     8  19   3   5  21   5   8   2  12
## ENSRNOG00000000008 Alx4 protein_coding    40  54  21  30  64  11  38  14  35
## ENSRNOG00000000009 Tmco5b protein_coding   0   0   0   0   2   0   0   0   0
## ENSRNOG00000000010 Cbln1 protein_coding    0   1   4   2   1   6   3   3   0
## ENSRNOG00000000012 Tcf15 protein_coding  482 501 322 319 906 598 834 630 619
##                                           26  27  28  29   3  30  31  32 33 34
## ENSRNOG00000000001 Arsj protein_coding    12  27  16  16  16  10  24  26 14 12
## ENSRNOG00000000007 Gad1 protein_coding     5  11  10   5  11   9   9   6 34 21
## ENSRNOG00000000008 Alx4 protein_coding    44  62  29  31  32  34  28  22 14 10
## ENSRNOG00000000009 Tmco5b protein_coding   0   0   0   0   0   0   0   0 18  7
## ENSRNOG00000000010 Cbln1 protein_coding    0   0   2   3   1   0   4   0 23 11
## ENSRNOG00000000012 Tcf15 protein_coding  648 867 758 472 493 296 468 289 19 11
##                                          35 36 37 38 39   4 40 41 42 43 44 45
## ENSRNOG00000000001 Arsj protein_coding   30 13 28 35 22  13 24 16 17 29 27 36
## ENSRNOG00000000007 Gad1 protein_coding   41 34 37 42 28   6 17 25 18 44 35 41
## ENSRNOG00000000008 Alx4 protein_coding   17  6 17 25 13  18 16  6  8 40 15 17
## ENSRNOG00000000009 Tmco5b protein_coding 20 11 23 29 26   0 22  5  8 20 10 22
## ENSRNOG00000000010 Cbln1 protein_coding  24 16 30 18 12   0 21 17 14 24 13 20
## ENSRNOG00000000012 Tcf15 protein_coding  10  6  5 18  2 247  7  4  6 19 14 13
##                                          46 47 48 49   5 50 51 52 53 54 55 56
## ENSRNOG00000000001 Arsj protein_coding   14  8 26 25  19 19 37 31 31 20 14 16
## ENSRNOG00000000007 Gad1 protein_coding   38 41 32 40   8 49 52 48 33 42 26 38
## ENSRNOG00000000008 Alx4 protein_coding    9 31  6 19  25 17 25 18 13 14 25  5
## ENSRNOG00000000009 Tmco5b protein_coding 13 21  7 22   0 27 36 33 22 33 24 17
## ENSRNOG00000000010 Cbln1 protein_coding  25 21 13 18   0 25 37 16 24 23 24 24
## ENSRNOG00000000012 Tcf15 protein_coding  10  6 11  4 832 15  8 15 10 15  6  7
##                                          57 58 59   6 60 61 62 63 64   7   8
## ENSRNOG00000000001 Arsj protein_coding   21 26 14  18 24 32 33 30 30  14  15
## ENSRNOG00000000007 Gad1 protein_coding   58 38 34  19 49 44 53 48 41   4   3
## ENSRNOG00000000008 Alx4 protein_coding   20 12 12  28 22 27 15 18 21  14  14
## ENSRNOG00000000009 Tmco5b protein_coding 28 15 15   0 10 23 21 25 18   0   0
## ENSRNOG00000000010 Cbln1 protein_coding  17 27 27   0 25 27 26 30 32   0   2
## ENSRNOG00000000012 Tcf15 protein_coding  12  8  5 646 10 10  9 10  4 351 269
##                                            9 95 96
## ENSRNOG00000000001 Arsj protein_coding    13  0  0
## ENSRNOG00000000007 Gad1 protein_coding     6  0  0
## ENSRNOG00000000008 Alx4 protein_coding    25  0  0
## ENSRNOG00000000009 Tmco5b protein_coding   0  0  0
## ENSRNOG00000000010 Cbln1 protein_coding    0  0  0
## ENSRNOG00000000012 Tcf15 protein_coding  438  3 33
write.table(x,file="countmatrix_skeletal.tsv",quote=FALSE,sep="\t")

y <- x/colSums(x)*1000000

write.table(y,file="rpmmatrix_skeletal.tsv",quote=FALSE,sep="\t")

Samplesheet

samplesheet <- read.table("samplesheet.tsv",header=TRUE)
samplesheet$UDI <- gsub("UDI_0","",samplesheet$UDI)
samplesheet$UDI <- gsub("UDI_","",samplesheet$UDI)
samplesheet$label <- paste(samplesheet$fraction,samplesheet$Group)
samplesheet <- samplesheet[order(samplesheet$UDI),]
samplesheet %>% kbl() %>% kable_paper("hover", full_width = F)
UDI fraction Group RatID label
1 1 wholetissue T 1 wholetissue T
10 10 wholetissue T 10 wholetissue T
11 11 wholetissue S 11 wholetissue S
12 12 wholetissue S 12 wholetissue S
13 13 wholetissue T 13 wholetissue T
14 14 wholetissue T 14 wholetissue T
15 15 wholetissue S 15 wholetissue S
16 16 wholetissue S 16 wholetissue S
17 17 wholetissue T 17 wholetissue T
18 18 wholetissue T 18 wholetissue T
19 19 wholetissue S 19 wholetissue S
2 2 wholetissue T 2 wholetissue T
20 20 wholetissue S 20 wholetissue S
21 21 wholetissue T 21 wholetissue T
22 22 wholetissue T 22 wholetissue T
23 23 wholetissue S 23 wholetissue S
24 24 wholetissue S 24 wholetissue S
25 25 wholetissue T 25 wholetissue T
26 26 wholetissue T 26 wholetissue T
27 27 wholetissue T 27 wholetissue T
28 28 wholetissue T 28 wholetissue T
29 29 wholetissue S 29 wholetissue S
3 3 wholetissue S 3 wholetissue S
30 30 wholetissue S 30 wholetissue S
31 31 wholetissue S 31 wholetissue S
32 32 wholetissue S 32 wholetissue S
33 33 mito T 1 mito T
34 34 mito T 2 mito T
35 35 mito S 3 mito S
36 36 mito S 4 mito S
37 37 mito T 5 mito T
38 38 mito T 6 mito T
39 39 mito S 7 mito S
4 4 wholetissue S 4 wholetissue S
40 40 mito S 8 mito S
41 41 mito T 9 mito T
42 42 mito T 10 mito T
43 43 mito S 11 mito S
44 44 mito S 12 mito S
45 45 mito T 13 mito T
46 46 mito T 14 mito T
47 47 mito S 15 mito S
48 48 mito S 16 mito S
49 49 mito T 17 mito T
5 5 wholetissue T 5 wholetissue T
50 50 mito T 18 mito T
51 51 mito S 19 mito S
52 52 mito S 20 mito S
53 53 mito T 21 mito T
54 54 mito T 22 mito T
55 55 mito S 23 mito S
56 56 mito S 24 mito S
57 57 mito T 25 mito T
58 58 mito T 26 mito T
59 59 mito T 27 mito T
6 6 wholetissue T 6 wholetissue T
60 60 mito T 28 mito T
61 61 mito S 29 mito S
62 62 mito S 30 mito S
63 63 mito S 31 mito S
64 64 mito S 32 mito S
7 7 wholetissue S 7 wholetissue S
8 8 wholetissue S 8 wholetissue S
9 9 wholetissue T 9 wholetissue T

Overall clustering with multidimensional scaling

Sample 95 is a HUMAN mito control sample. Sample 96 is a RAT mito control sample. These should be excluded from main results.

ss <- samplesheet
ss <- ss[order(ss$UDI),]
colours = c('pink', 'lightblue','lightgreen','gray')
mds <- cmdscale(dist(t(x)))
XMAX=max(mds[,1])*1.1
XMIN=min(mds[,1])*1.1
plot( mds*1.05 , cex=2, pch=19, xlab="Coordinate 1", ylab="Coordinate 2",
  col = colours[as.factor(ss$label)] ,  type = "p" ,
  xlim=c(XMIN,XMAX),main="MDS plot",bty="n")
text(mds, labels=colnames(x) )
legend('topright', col=colours, legend=levels(as.factor(ss$label)), pch = 16, cex = 1)

x <- x[,which(! colnames(x) %in% c("95","96"))]

colours = c('pink', 'lightblue','lightgreen','gray')
mds <- cmdscale(dist(t(x)))
XMAX=max(mds[,1])*1.1
XMIN=min(mds[,1])*1.1
plot( mds*1.05 , cex=2, pch=19, xlab="Coordinate 1", ylab="Coordinate 2",
  col = colours[as.factor(ss$label)] ,  type = "p" ,
  xlim=c(XMIN,XMAX),main="MDS plot",bty="n")
text(mds, labels=colnames(x) )
legend('topright', col=colours, legend=levels(as.factor(ss$label)), pch = 16, cex = 1)

ss$nreads<-colSums(x)
par(mar=c(5,10,5,3))
barplot(colSums(x),horiz=TRUE,las=2,main="number of reads per sample",cex.names=0.5)

par(mai=c(1.02,0.82,0.82,0.42))

Check purity of mito fraction samples

Here I’m quantifying the mitochondrial read fraction. That is the number of mt reads divided by the total number of reads. We can see that purity is highly variable.

par(mar=c(5,10,5,3))

mtgenes <- c("ENSRNOG00000029042",
"ENSRNOG00000029070",
"ENSRNOG00000029145",
"ENSRNOG00000029171",
"ENSRNOG00000029301",
"ENSRNOG00000029389",
"ENSRNOG00000029677",
"ENSRNOG00000029707",
"ENSRNOG00000029954",
"ENSRNOG00000029971",
"ENSRNOG00000030339",
"ENSRNOG00000030371",
"ENSRNOG00000030478",
"ENSRNOG00000030644",
"ENSRNOG00000030700",
"ENSRNOG00000031033",
"ENSRNOG00000031053",
"ENSRNOG00000031333",
"ENSRNOG00000031667",
"ENSRNOG00000031685",
"ENSRNOG00000031766",
"ENSRNOG00000031780",
"ENSRNOG00000031979",
"ENSRNOG00000032112",
"ENSRNOG00000032274",
"ENSRNOG00000032320",
"ENSRNOG00000032578",
"ENSRNOG00000032609",
"ENSRNOG00000032882",
"ENSRNOG00000032997",
"ENSRNOG00000033299",
"ENSRNOG00000033545",
"ENSRNOG00000033615",
"ENSRNOG00000033932",
"ENSRNOG00000033957",
"ENSRNOG00000034234",
"ENSRNOG00000043866")

mtcounts <-  lapply(mtgenes, function(i) { x[grep(i,rownames(x)),] } )
mtcounts <- do.call(rbind,mtcounts)
mtfrac<- colSums(mtcounts)/colSums(x)
barplot(mtfrac,horiz=TRUE,las=2,main="Proportion mitochondrial reads",cex.names=1)

par(mai=c(1.02,0.82,0.82,0.42))

mylevels <- levels(as.factor(ss$label))
mylevels
## [1] "mito S"        "mito T"        "wholetissue S" "wholetissue T"
y <- x[,which(ss$label==mylevels[1])]
mitoS <- colSums(y[grep("^MT",rownames(y)),]) / colSums(y)
y <- x[,which(ss$label==mylevels[2])] 
mitoT <- colSums(y[grep("^MT",rownames(y)),]) / colSums(y)
y <- x[,which(ss$label==mylevels[3])] 
wholeS <- colSums(y[grep("^MT",rownames(y)),]) / colSums(y)
y <- x[,which(ss$label==mylevels[4])] 
wholeT <- colSums(y[grep("^MT",rownames(y)),]) / colSums(y)

boxplot(mitoS,mitoT,wholeS,wholeT,names=mylevels,ylab="mito frac")

median(wholeS)
## [1] 0
median(wholeT)
## [1] 0
table(mitoS>0.4)
## 
## FALSE 
##    16
table(mitoT>0.4)
## 
## FALSE 
##    16
ss$mtfrac <- mtfrac

ss %>% kbl() %>% kable_paper("hover", full_width = F)
UDI fraction Group RatID label nreads mtfrac
1 1 wholetissue T 1 wholetissue T 22353545 0.2350623
10 10 wholetissue T 10 wholetissue T 32594141 0.3531781
11 11 wholetissue S 11 wholetissue S 28160771 0.3290185
12 12 wholetissue S 12 wholetissue S 30053603 0.2640899
13 13 wholetissue T 13 wholetissue T 29989340 0.2972082
14 14 wholetissue T 14 wholetissue T 27775661 0.3282823
15 15 wholetissue S 15 wholetissue S 25070009 0.3071025
16 16 wholetissue S 16 wholetissue S 27766116 0.3024493
17 17 wholetissue T 17 wholetissue T 30401193 0.3348082
18 18 wholetissue T 18 wholetissue T 33702596 0.3547981
19 19 wholetissue S 19 wholetissue S 30788467 0.3233271
2 2 wholetissue T 2 wholetissue T 19477976 0.3181207
20 20 wholetissue S 20 wholetissue S 30589099 0.3209945
21 21 wholetissue T 21 wholetissue T 39401323 0.3400537
22 22 wholetissue T 22 wholetissue T 27353327 0.3221241
23 23 wholetissue S 23 wholetissue S 32055355 0.3547227
24 24 wholetissue S 24 wholetissue S 28630316 0.2413714
25 25 wholetissue T 25 wholetissue T 34339726 0.2938285
26 26 wholetissue T 26 wholetissue T 28163757 0.2587561
27 27 wholetissue T 27 wholetissue T 30543678 0.3380834
28 28 wholetissue T 28 wholetissue T 30906624 0.2776614
29 29 wholetissue S 29 wholetissue S 24309444 0.3591212
3 3 wholetissue S 3 wholetissue S 23906716 0.3067514
30 30 wholetissue S 30 wholetissue S 26615136 0.2724498
31 31 wholetissue S 31 wholetissue S 27019886 0.3536612
32 32 wholetissue S 32 wholetissue S 20876971 0.3177831
33 33 mito T 1 mito T 5474018 0.8850587
34 34 mito T 2 mito T 9988027 0.9542655
35 35 mito S 3 mito S 1469458 0.4723449
36 36 mito S 4 mito S 1168009 0.5294668
37 37 mito T 5 mito T 1353729 0.4835828
38 38 mito T 6 mito T 914913 0.1714841
39 39 mito S 7 mito S 1780282 0.6982787
4 4 wholetissue S 4 wholetissue S 25823475 0.2520651
40 40 mito S 8 mito S 712067 0.2883577
41 41 mito T 9 mito T 581457 0.1897888
42 42 mito T 10 mito T 530397 0.1889811
43 43 mito S 11 mito S 1700043 0.4096285
44 44 mito S 12 mito S 720111 0.2030354
45 45 mito T 13 mito T 971953 0.2654881
46 46 mito T 14 mito T 1086442 0.4345745
47 47 mito S 15 mito S 8344505 0.9222841
48 48 mito S 16 mito S 1296403 0.5776151
49 49 mito T 17 mito T 1010393 0.2873624
5 5 wholetissue T 5 wholetissue T 30755504 0.3586038
50 50 mito T 18 mito T 1316551 0.3568658
51 51 mito S 19 mito S 1483811 0.2887012
52 52 mito S 20 mito S 3085085 0.7437176
53 53 mito T 21 mito T 1628001 0.5591495
54 54 mito T 22 mito T 1252520 0.3879571
55 55 mito S 23 mito S 888159 0.1293935
56 56 mito S 24 mito S 2133695 0.6219164
57 57 mito T 25 mito T 1841975 0.5667390
58 58 mito T 26 mito T 3802593 0.7906639
59 59 mito T 27 mito T 7532868 0.9057456
6 6 wholetissue T 6 wholetissue T 26148398 0.3245408
60 60 mito T 28 mito T 6001091 0.8726930
61 61 mito S 29 mito S 1785865 0.5147013
62 62 mito S 30 mito S 1474773 0.4122207
63 63 mito S 31 mito S 1503193 0.4336855
64 64 mito S 32 mito S 4091353 0.8123151
7 7 wholetissue S 7 wholetissue S 25564135 0.2969021
8 8 wholetissue S 8 wholetissue S 22385138 0.2206553
9 9 wholetissue T 9 wholetissue T 24737507 0.3535222 
mito <- subset(ss,fraction=="mito")

plot(mito$nreads,mito$mtfrac,pch=19,col="lightblue",xlab="No. reads", ylab="mtDNA proportion",cex=3)
text(mito$nreads,mito$mtfrac,labels=mito$RatID)
abline(h=0.4,v=6E6,col="red")

plot(mito$nreads,mito$mtfrac,pch=19,col="lightblue",xlab="No. reads", ylab="mtDNA proportion",cex=3)
text(mito$nreads,mito$mtfrac,labels=mito$Group)
abline(h=0.4,v=6E6,col="red")

Functions

run_de <- function(ss,xx){

y <- round(xx)

# MDS
colours = c('yellow', 'orange')
mds <- cmdscale(dist(t(y)))
XMAX=max(mds[,1])*1.1
XMIN=min(mds[,1])*1.1
plot( mds*1.05 , cex=2 , pch=19, xlab="Coordinate 1", ylab="Coordinate 2",
  col = colours[as.factor(ss$trt)] ,  type = "p" , 
  xlim=c(XMIN,XMAX),main="MDS plot",bty="n")
text(mds, labels=colnames(y) )
legend('topright', col=colours, legend=c("ctrl","trt"), pch = 16, cex = 1.5)

# DE
dds <- DESeqDataSetFromMatrix(countData=y, colData = ss, design = ~ trt)
dds <- DESeq(dds)
de <- DESeq2::results(dds)
de <- de[order(de$pvalue),]
up <- rownames(subset(de, log2FoldChange>0 & padj<0.05 ))
dn <- rownames(subset(de, log2FoldChange<0 & padj<0.05 ))
str(up)
str(dn)

# MA plot
sig <-subset(de, padj < 0.05 )
GENESUP <- length(up)
GENESDN <- length(dn)
SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down")
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="smear plot")
points(log2(sig$baseMean),sig$log2FoldChange,
       pch=19, cex=0.5, col="red")
mtext(SUBHEADER)

# heatmap
yn <- y/colSums(y)*1000000
yf <- yn[which(rownames(yn) %in% rownames(de)[1:50]),]
mycols <- gsub("0","yellow",ss$trt)
mycols <- gsub("1","orange",mycols)
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(  as.matrix(yf), col=colfunc(25),scale="row",
    ColSideColors =mycols ,trace="none",
    margin = c(10,15), cexRow=0.6, cexCol=0.8 , main="Top 50 genes by p-val")
mtext("yellow=ctrl, orange=trt")

de <- merge(as.data.frame(de),yn,by=0)
rownames(de) <- de[,1]
de[,1]=NULL
return(de)
}

# covariate
run_de_cov <- function(ss,xx){

y <- round(xx)

# MDS
colours = c('yellow', 'orange')
mds <- cmdscale(dist(t(y)))
XMAX=max(mds[,1])*1.1
XMIN=min(mds[,1])*1.1
plot( mds*1.05 , cex=2 , pch=19, xlab="Coordinate 1", ylab="Coordinate 2",
  col = colours[as.factor(ss$trt)] ,  type = "p" ,
  xlim=c(XMIN,XMAX),main="MDS plot",bty="n")
text(mds, labels=colnames(y) )
legend('topright', col=colours, legend=c("ctrl","trt"), pch = 16, cex = 1.5)

# DE
dds <- DESeqDataSetFromMatrix(countData=y, colData = ss, design = ~ mtfrac + trt)
dds <- DESeq(dds)
de <- DESeq2::results(dds)
de <- de[order(de$pvalue),]
up <- rownames(subset(de, log2FoldChange>0 & padj<0.05 ))
dn <- rownames(subset(de, log2FoldChange<0 & padj<0.05 ))
str(up)
str(dn)

# MA plot
sig <-subset(de, padj < 0.05 )
GENESUP <- length(up)
GENESDN <- length(dn)
SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down")
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="smear plot")
points(log2(sig$baseMean),sig$log2FoldChange,
       pch=19, cex=0.5, col="red")
mtext(SUBHEADER)

# heatmap
yn <- y/colSums(y)*1000000
yf <- yn[which(rownames(yn) %in% rownames(de)[1:50]),]
mycols <- gsub("0","yellow",ss$trt)
mycols <- gsub("1","orange",mycols)
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(  as.matrix(yf), col=colfunc(25),scale="row",
    ColSideColors =mycols ,trace="none",
    margin = c(10,15), cexRow=0.6, cexCol=0.8 , main="Top 50 genes by p-val")
mtext("yellow=ctrl, orange=trt")

de <- merge(as.data.frame(de),yn,by=0)
rownames(de) <- de[,1]
de[,1]=NULL
return(de)
}

Set up contrasts

  1. Whole tissue S versus T

  2. Mito fraction S versus T

  3. Whole versus mito fraction in S

  4. Whole versus mito fraction in T

ss1 <- subset(samplesheet,fraction=="wholetissue")
ss1$trt <- grepl("T",ss1$Group)*1
x1 <- x[,which(colnames(x) %in% rownames(ss1))]

ss2 <- subset(samplesheet,fraction=="mito")
ss2$trt <- grepl("T",ss2$Group)*1
x2 <- x[,which(colnames(x) %in% rownames(ss2))]

ss3 <- subset(samplesheet,Group=="S")
ss3$trt <- grepl("mito",ss3$fraction)*1
x3 <- x[,which(colnames(x) %in% rownames(ss3))]

ss4 <- subset(samplesheet,Group=="T")
ss4$trt <- grepl("mito",ss4$fraction)*1
x4 <- x[,which(colnames(x) %in% rownames(ss4))]

ss5 <- subset(mito,nreads>6E6 & mtfrac>0.4)
ss5$trt <- grepl("T",ss5$Group)*1
x5 <- x[,which(colnames(x) %in% rownames(ss5))]

DE

Here, were using DESeq2 to perform differential expression analysis for the specified contrasts. The run_de function does the analysis and generate the charts. Here we actually run the analysis.

  1. Whole tissue S versus T: 1 DEG

  2. Mito fraction S versus T: 0 DEGs

  3. Whole versus mito fraction in S: 22754 DEGs

  4. Whole versus mito fraction in T: 22406 DEGs

These results suggest that:

  • The differences between S and T are very subtle.

  • Intra-group variation is fairly large.

ss1 %>% kbl() %>% kable_paper("hover", full_width = F)
UDI fraction Group RatID label trt
1 1 wholetissue T 1 wholetissue T 1
10 10 wholetissue T 10 wholetissue T 1
11 11 wholetissue S 11 wholetissue S 0
12 12 wholetissue S 12 wholetissue S 0
13 13 wholetissue T 13 wholetissue T 1
14 14 wholetissue T 14 wholetissue T 1
15 15 wholetissue S 15 wholetissue S 0
16 16 wholetissue S 16 wholetissue S 0
17 17 wholetissue T 17 wholetissue T 1
18 18 wholetissue T 18 wholetissue T 1
19 19 wholetissue S 19 wholetissue S 0
2 2 wholetissue T 2 wholetissue T 1
20 20 wholetissue S 20 wholetissue S 0
21 21 wholetissue T 21 wholetissue T 1
22 22 wholetissue T 22 wholetissue T 1
23 23 wholetissue S 23 wholetissue S 0
24 24 wholetissue S 24 wholetissue S 0
25 25 wholetissue T 25 wholetissue T 1
26 26 wholetissue T 26 wholetissue T 1
27 27 wholetissue T 27 wholetissue T 1
28 28 wholetissue T 28 wholetissue T 1
29 29 wholetissue S 29 wholetissue S 0
3 3 wholetissue S 3 wholetissue S 0
30 30 wholetissue S 30 wholetissue S 0
31 31 wholetissue S 31 wholetissue S 0
32 32 wholetissue S 32 wholetissue S 0
4 4 wholetissue S 4 wholetissue S 0
5 5 wholetissue T 5 wholetissue T 1
6 6 wholetissue T 6 wholetissue T 1
7 7 wholetissue S 7 wholetissue S 0
8 8 wholetissue S 8 wholetissue S 0
9 9 wholetissue T 9 wholetissue T 1 
de1 <- run_de(ss1,x1)
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 214 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
##  chr [1:15] "ENSRNOG00000019961 Tcte1 protein_coding" ...
##  chr "ENSRNOG00000016581 Serpinb1a protein_coding"

as.data.frame(de1[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)
baseMean log2FoldChange lfcSE stat pvalue padj 1 10 11 12 13 14 15 16 17 18 19 2 20 21 22 23 24 25 26 27 28 29 3 30 31 32 4 5 6 7 8 9
ENSRNOG00000000001 Arsj protein_coding 19.5250941 0.1851575 0.1862240 0.9942731 0.3200899 0.7983311 0.8052414 1.4315403 0.6710345 0.7157701 1.3868046 0.8052414 0.5368276 0.6710345 2.1920460 0.8947127 0.8947127 0.9841839 0.7157701 0.7157701 0.8052414 1.0736552 1.2973334 0.8947127 0.5368276 1.2078621 0.7157701 0.7157701 0.7157701 0.4473563 1.0736552 1.1631265 0.5815632 0.8499770 0.8052414 0.6262989 0.6710345 0.5815632
ENSRNOG00000000007 Gad1 protein_coding 7.6926925 0.1888669 0.2841229 0.6647367 0.5062189 0.8818083 0.2454429 0.1534018 0.3068036 0.1840822 0.3681643 0.0306804 0.1840822 0.0920411 0.2761232 0.2454429 0.5829269 0.0920411 0.1534018 0.6442876 0.1534018 0.2454429 0.0613607 0.3681643 0.1534018 0.3374840 0.3068036 0.1534018 0.3374840 0.2761232 0.2761232 0.1840822 0.1840822 0.2454429 0.5829269 0.1227214 0.0920411 0.1840822
ENSRNOG00000000008 Alx4 protein_coding 31.2722329 0.1822232 0.2237965 0.8142362 0.4155096 0.8440259 0.2130623 1.6334780 1.6334780 0.5681663 0.8877598 1.5624572 1.3849053 0.8167390 2.5922586 1.4204157 1.9175611 0.7457182 1.0653117 2.2726650 0.3906143 1.3493949 0.4971455 1.2428637 1.5624572 2.2016443 1.0298014 1.1008221 1.1363325 1.2073533 0.9942910 0.7812286 0.6391870 0.8877598 0.9942910 0.4971455 0.4971455 0.8877598
ENSRNOG00000000009 Tmco5b protein_coding 0.1412217 0.7722730 2.0757455 0.3720461 0.7098585 NA 0.0332739 0.0332739 0.0000000 0.0000000 0.0000000 0.0332739 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0665478 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
ENSRNOG00000000010 Cbln1 protein_coding 1.1882585 -0.5604697 0.7603507 -0.7371199 0.4610494 0.8623395 0.0000000 0.0000000 0.0333452 0.0333452 0.0000000 0.0333452 0.0000000 0.0333452 0.0333452 0.0000000 0.0333452 0.1333807 0.0666904 0.0333452 0.2000711 0.1000355 0.1000355 0.0000000 0.0000000 0.0000000 0.0666904 0.1000355 0.0333452 0.0000000 0.1333807 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0666904 0.0000000
ENSRNOG00000000012 Tcf15 protein_coding 488.8633289 0.3824149 0.1242079 3.0788299 0.0020782 0.2114262 10.0087627 24.4098601 12.0969218 18.9734459 25.2739260 16.0572236 11.0528423 12.4929520 13.2490096 17.3533224 18.0373745 11.5928834 11.4848752 32.6184857 21.5296406 30.0262881 22.6817284 22.2856983 23.3297778 31.2143787 27.2900796 16.9932950 17.7493526 10.6568121 16.8492840 10.4047929 8.8926777 29.9542826 23.2577723 12.6369630 9.6847380 15.7692017
ENSRNOG00000000017 Steap1 protein_coding 0.6812740 -1.6951967 1.0207209 -1.6607838 0.0967569 NA 0.0000000 0.0000000 0.0797766 0.0797766 0.0398883 0.0000000 0.0000000 0.0000000 0.1196649 0.0000000 0.0398883 0.0398883 0.0000000 0.0000000 0.0000000 0.0797766 0.1595532 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1196649 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1196649 0.0000000
ENSRNOG00000000021 AABR07061902.1 protein_coding 0.0286358 0.0509230 2.9741591 0.0171218 0.9863394 NA 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0360151 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
ENSRNOG00000000024 Hebp1 protein_coding 59.0108585 0.0369259 0.0997940 0.3700210 0.7113668 0.9420024 1.3157378 3.0261970 2.5656888 1.9407133 2.7959429 1.7104592 1.6117789 1.7104592 2.4999019 2.0065002 1.8749264 1.1841641 2.4341150 2.4999019 2.4670084 2.6314757 2.4341150 1.8749264 2.1380740 2.5656888 2.4012216 1.6775657 1.7762461 1.5788854 2.0065002 1.0525903 1.2828444 1.8091395 2.2367543 2.2696478 1.2170575 1.2828444
ENSRNOG00000000033 Tmcc2 protein_coding 1002.6860405 0.0004456 0.0798111 0.0055826 0.9955458 0.9989304 31.4219118 29.8790040 23.0842752 35.4275380 33.9143014 29.3745918 26.3481187 29.0185361 31.3922405 35.3978667 30.0570318 18.9302925 34.1220006 36.3770197 32.1933658 24.8052109 39.7595485 40.2342894 36.1989919 39.7298772 39.1364511 25.0129100 27.5052996 33.2615327 23.3809882 14.7466385 32.9054771 27.8910266 27.8020126 31.8669814 29.9086753 18.9896351
ENSRNOG00000000034 Nuak2 protein_coding 59.9809533 0.1199926 0.1371948 0.8746146 0.3817836 0.8278052 1.9812614 2.6958146 2.3385380 2.1436598 2.7607740 1.7214238 1.7863832 1.6239847 2.0137411 1.2342284 1.7539035 1.2667081 1.3641472 1.9487817 2.4034974 3.3454085 3.1505304 2.8257334 2.2735786 3.7026852 1.3316675 1.5265456 1.3316675 1.0718299 2.1436598 1.2342284 1.4615862 2.9881319 1.7539035 2.1436598 1.0393502 1.7539035
ENSRNOG00000000035 Tmed11 protein_coding 0.0000000 NA NA NA NA NA 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
ENSRNOG00000000036 Klhdc8a protein_coding 2.5911505 -0.2780132 0.5489821 -0.5064157 0.6125648 0.9165793 0.0326914 0.2288397 0.0000000 0.1307655 0.0326914 0.1961483 0.0653828 0.0000000 0.2615311 0.0000000 0.0653828 0.0000000 0.0000000 0.0653828 0.0980742 0.0653828 0.0326914 0.1634569 0.0653828 0.0653828 0.0326914 0.0980742 0.1961483 0.0980742 0.1307655 0.0326914 0.0326914 0.0000000 0.0000000 0.0653828 0.3596052 0.0653828
ENSRNOG00000000040 RGD1304622 protein_coding 126.7630749 -0.0264477 0.0800284 -0.3304793 0.7410379 0.9488852 2.9186837 4.2891961 3.0709629 4.1876766 3.9592579 3.0963427 1.9288692 3.7308392 3.1978622 3.7815989 2.8425441 2.1065283 3.7562190 3.7815989 2.9186837 2.7410247 4.5429947 3.1724823 3.4262809 3.7815989 3.4262809 2.6648851 2.9694434 3.2486219 3.0202032 1.8781095 3.5278003 4.1876766 3.6546996 3.3755212 3.0202032 2.5633657
ENSRNOG00000000041 Slc26a1 protein_coding 1.2304535 2.0624148 1.0491798 1.9657400 0.0493286 0.5348888 0.4021449 0.0365586 0.0000000 0.0000000 0.0000000 0.0365586 0.0000000 0.0365586 0.0000000 0.0365586 0.0731172 0.2924690 0.0000000 0.0000000 0.1096759 0.0000000 0.0000000 0.0000000 0.0000000 0.0365586 0.0000000 0.0000000 0.1096759 0.0365586 0.0000000 0.0000000 0.0000000 0.0000000 0.0365586 1.0602001 0.0000000 0.0000000
ENSRNOG00000000042 Xpr1 protein_coding 469.5901686 -0.0828081 0.0473421 -1.7491414 0.0802666 0.5971745 14.4437645 16.0035663 16.2843307 16.5962910 15.7228020 14.0694121 12.9151588 14.5997447 15.7228020 16.3155267 17.3449959 10.2011037 15.1612734 20.6517757 15.7539980 15.2236654 18.3120730 17.5321721 15.9723703 15.8787822 17.1266236 12.1352579 13.1335310 15.0676853 12.2912381 13.1647271 14.5373527 12.7903747 13.1335310 13.9446280 12.3536302 10.9498085
ENSRNOG00000000043 Idua protein_coding 357.0035160 -0.0150644 0.0628613 -0.2396445 0.8106058 0.9641349 10.4085474 13.1678602 13.8664205 14.4601967 15.1238289 11.7358118 12.5391560 13.8314925 12.0501639 14.0759885 11.1769636 9.9544832 10.5133314 14.8094768 13.3425003 11.6659558 12.9582922 15.6477491 12.7836521 16.2764533 15.6128210 10.8975395 9.1860670 12.6788681 13.1329322 10.6181154 11.0372516 13.9362765 12.1549479 11.2118916 12.8185801 9.6401311
ENSRNOG00000000044 Tmem175 protein_coding 263.6431323 0.0553128 0.0570255 0.9699666 0.3320631 0.8053578 6.6104197 9.3477741 8.7944790 7.7752513 8.2994256 7.8043721 6.9016276 7.6587682 7.8334929 8.1538216 9.2021701 5.0378969 6.2900910 11.9977661 7.7752513 9.4933780 8.9983246 11.1823839 7.3384395 10.2505186 9.8428275 6.4648157 5.5620712 7.7752513 6.7269028 5.9697623 6.9889900 7.8043721 7.5714058 6.2900910 6.0571246 7.1345939
ENSRNOG00000000047 Cd82 protein_coding 143.2116013 0.0165392 0.0695378 0.2378447 0.8120016 0.9644264 4.4383283 5.7875801 5.5035271 6.7817657 5.2904873 4.4028217 4.2607952 4.0832620 7.1723385 5.3615006 4.4738349 3.7992090 5.0774476 8.1310175 5.2549807 7.1013253 6.2491663 5.2904873 4.8644078 6.4266994 5.8230867 4.2252886 5.6100470 4.3318084 4.6158614 3.5506626 4.1897819 6.6397392 4.6513681 5.0064343 3.7637024 4.1187687
ENSRNOG00000000048 Gak protein_coding 638.5153786 -0.0480836 0.0463336 -1.0377698 0.2993773 0.7896194 16.6646597 23.2453996 22.3941596 23.1799196 24.9151396 17.8432997 15.9443797 20.5607196 20.5607196 27.0432395 21.1827796 12.3102398 21.4119596 26.9450195 24.3912996 25.3080196 30.9392995 26.8140595 23.3108796 24.5549996 25.6354195 17.0575397 17.5486397 18.5963197 22.1322396 16.0753397 17.6795997 21.8375796 20.1678396 21.5429196 19.5785197 16.4027397 
write.table(de1,file="de1.tsv",quote=FALSE,sep="\t")

ss2 %>% kbl() %>% kable_paper("hover", full_width = F)
UDI fraction Group RatID label trt
33 33 mito T 1 mito T 1
34 34 mito T 2 mito T 1
35 35 mito S 3 mito S 0
36 36 mito S 4 mito S 0
37 37 mito T 5 mito T 1
38 38 mito T 6 mito T 1
39 39 mito S 7 mito S 0
40 40 mito S 8 mito S 0
41 41 mito T 9 mito T 1
42 42 mito T 10 mito T 1
43 43 mito S 11 mito S 0
44 44 mito S 12 mito S 0
45 45 mito T 13 mito T 1
46 46 mito T 14 mito T 1
47 47 mito S 15 mito S 0
48 48 mito S 16 mito S 0
49 49 mito T 17 mito T 1
50 50 mito T 18 mito T 1
51 51 mito S 19 mito S 0
52 52 mito S 20 mito S 0
53 53 mito T 21 mito T 1
54 54 mito T 22 mito T 1
55 55 mito S 23 mito S 0
56 56 mito S 24 mito S 0
57 57 mito T 25 mito T 1
58 58 mito T 26 mito T 1
59 59 mito T 27 mito T 1
60 60 mito T 28 mito T 1
61 61 mito S 29 mito S 0
62 62 mito S 30 mito S 0
63 63 mito S 31 mito S 0
64 64 mito S 32 mito S 0 
de2 <- run_de(ss2,x2)
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
## 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
##  chr(0) 
##  chr(0)

as.data.frame(de2[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)
baseMean log2FoldChange lfcSE stat pvalue padj 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
ENSRNOG00000000001 Arsj protein_coding 22.672617 -0.0672219 0.1464715 -0.4589421 0.6462757 0.9992526 2.557536 2.1921740 5.480435 2.374855 5.115073 6.393841 4.018986 4.384348 2.922899 3.105580 5.297754 4.932391 6.576522 2.557536 1.461449 4.749710 4.567029 3.470942 6.759203 5.663116 5.663116 3.653623 2.557536 2.922899 3.836304 4.749710 2.557536 4.384348 5.845797 6.028479 5.480435 5.480435
ENSRNOG00000000007 Gad1 protein_coding 36.450485 0.0696246 0.0996736 0.6985256 0.4848485 0.9992526 3.404076 2.1025173 4.104915 3.404076 3.704435 4.205035 2.803357 1.702038 2.502997 1.802158 4.405274 3.504196 4.104915 3.804555 4.104915 3.203836 4.004795 4.905874 5.206233 4.805754 3.303956 4.205035 2.603117 3.804555 5.806953 3.804555 3.404076 4.905874 4.405274 5.306353 4.805754 4.104915
ENSRNOG00000000008 Alx4 protein_coding 15.715072 -0.2137098 0.1801934 -1.1860019 0.2356216 0.9992526 9.527322 6.8052302 11.568891 4.083138 11.568891 17.013076 8.846799 10.888368 4.083138 5.444184 27.220921 10.207845 11.568891 6.124707 21.096214 4.083138 12.929937 11.568891 17.013076 12.249414 8.846799 9.527322 17.013076 3.402615 13.610461 8.166276 8.166276 14.971507 18.374122 10.207845 12.249414 14.290983
ENSRNOG00000000009 Tmco5b protein_coding 18.583695 -0.0666880 0.1629263 -0.4093139 0.6823093 0.9992526 15.410840 5.9931045 17.123156 9.417736 19.691629 24.828576 22.260102 18.835471 4.280789 6.849262 17.123156 8.561578 18.835471 11.130051 17.979313 5.993105 18.835471 23.116260 30.821680 28.253207 18.835471 28.253207 20.547787 14.554682 23.972418 12.842367 12.842367 8.561578 19.691629 17.979313 21.403945 15.410840
ENSRNOG00000000010 Cbln1 protein_coding 21.097175 0.0720593 0.1289348 0.5588821 0.5762422 0.9992526 16.990107 8.1257032 17.728807 11.819205 22.161009 13.296605 8.864403 15.512706 12.557905 10.341804 17.728807 9.603104 14.774006 18.467507 15.512706 9.603104 13.296605 18.467507 27.331911 11.819205 17.728807 16.990107 17.728807 17.728807 12.557905 19.944908 19.944908 18.467507 19.944908 19.206207 22.161009 23.638409
ENSRNOG00000000012 Tcf15 protein_coding 9.386898 0.3250874 0.2209877 1.4710658 0.1412733 0.9992526 20.767002 12.0230011 10.930001 6.558001 5.465001 19.674002 2.186000 7.651001 4.372000 6.558001 20.767002 15.302001 14.209001 10.930001 6.558001 12.023001 4.372000 16.395001 8.744001 16.395001 10.930001 16.395001 6.558001 7.651001 13.116001 8.744001 5.465001 10.930001 10.930001 9.837001 10.930001 4.372000
ENSRNOG00000000017 Steap1 protein_coding 9.063332 -0.2294720 0.2102640 -1.0913522 0.2751180 0.9992526 6.178796 2.2468351 8.425631 5.617088 11.234175 3.370253 2.246835 5.055379 3.931961 3.931961 9.549049 5.617088 6.740505 3.370253 4.493670 4.493670 6.740505 3.370253 7.302214 3.931961 3.931961 2.246835 6.178796 7.863923 6.178796 1.685126 3.370253 5.055379 5.617088 3.370253 6.740505 8.425631
ENSRNOG00000000021 AABR07061902.1 protein_coding 7.566953 0.1781270 0.2546045 0.6996223 0.4841632 0.9992526 4.213087 5.6174489 21.065433 7.021811 21.065433 11.234898 14.043622 14.043622 11.234898 5.617449 23.874158 9.830536 16.852347 11.234898 4.213087 8.426173 12.639260 9.830536 16.852347 8.426173 4.213087 12.639260 8.426173 5.617449 18.256709 21.065433 11.234898 2.808724 15.447985 8.426173 7.021811 0.000000
ENSRNOG00000000024 Hebp1 protein_coding 9.559144 -0.2825425 0.1915226 -1.4752432 0.1401472 0.9992526 8.599088 10.3189058 15.478359 6.879270 10.318906 12.038723 12.038723 12.038723 15.478359 5.159453 39.555805 18.917994 22.357629 15.478359 12.038723 15.478359 8.599088 18.917994 17.198176 18.917994 10.318906 27.517082 24.077447 25.797265 18.917994 18.917994 17.198176 18.917994 25.797265 25.797265 30.956717 18.917994
ENSRNOG00000000033 Tmcc2 protein_coding 35.716146 0.0056460 0.1023698 0.0551533 0.9560163 0.9992526 58.446786 43.3637445 69.759067 67.873687 90.498249 77.300588 49.019885 65.988307 35.822224 60.332166 98.039770 65.988307 71.644447 73.529828 52.790645 49.019885 62.217546 86.727489 115.008192 71.644447 73.529828 75.415208 65.988307 60.332166 81.071348 79.185968 58.446786 41.478364 86.727489 67.873687 103.695911 77.300588
ENSRNOG00000000034 Nuak2 protein_coding 33.921616 0.0656914 0.1042057 0.6304015 0.5284319 0.9992526 19.999494 14.1172900 18.234833 13.529070 15.881951 31.763902 14.705510 10.587967 12.940849 9.999747 32.352123 18.234833 20.587715 19.411274 17.058392 12.940849 19.411274 26.469919 33.528564 29.411021 14.705510 27.058139 24.705257 21.175935 23.528817 27.646359 23.528817 22.352376 24.705257 31.763902 20.587715 22.352376
ENSRNOG00000000035 Tmed11 protein_coding 4.302859 0.2769968 0.2775349 0.9980611 0.3182497 0.9992526 6.943374 8.3320488 6.943374 6.943374 4.166024 2.777350 5.554699 5.554699 2.777350 5.554699 5.554699 1.388675 8.332049 4.166024 4.166024 2.777350 8.332049 13.886748 9.720724 6.943374 5.554699 6.943374 6.943374 5.554699 9.720724 2.777350 5.554699 6.943374 6.943374 5.554699 5.554699 6.943374
ENSRNOG00000000036 Klhdc8a protein_coding 30.737847 -0.0251103 0.1196827 -0.2098075 0.8338179 0.9992526 31.894546 14.4039887 23.663696 28.807977 23.663696 45.269679 26.750265 14.403989 23.663696 13.375132 79.221938 34.981115 30.865690 26.750265 29.836834 18.519414 40.125397 36.009972 56.587098 31.894546 45.269679 34.981115 33.952259 37.038828 42.183110 38.067684 32.923403 28.807977 38.067684 43.211966 44.240822 28.807977
ENSRNOG00000000040 RGD1304622 protein_coding 50.075230 -0.0891120 0.1078258 -0.8264443 0.4085521 0.9992526 55.226142 20.2495853 47.862656 48.783092 41.419606 44.180914 27.613071 37.737864 34.056121 33.135685 86.520956 38.658299 52.464835 35.896992 55.226142 48.783092 49.703528 43.260478 84.680084 33.135685 23.931328 56.146577 49.703528 37.737864 68.112242 49.703528 57.067013 46.021785 61.669192 57.067013 59.828320 54.305706
ENSRNOG00000000041 Slc26a1 protein_coding 30.023635 -0.0378155 0.1162928 -0.3251753 0.7450484 0.9992526 3.834859 0.8388754 5.632449 2.636466 3.715020 4.194377 2.636466 2.876144 2.756305 2.636466 4.793574 2.516626 4.673734 3.235662 4.553895 3.235662 2.756305 5.153092 7.669718 3.235662 4.194377 2.995984 4.074538 2.756305 4.434056 5.153092 2.995984 3.595180 4.793574 5.153092 4.074538 3.595180
ENSRNOG00000000042 Xpr1 protein_coding 37.987794 0.1199608 0.1148738 1.0442837 0.2963542 0.9992526 17.741397 19.2841269 29.311873 28.540508 33.168698 33.168698 20.055492 17.741397 21.598222 16.970032 40.110984 26.226413 26.997778 25.455047 21.598222 18.512762 33.940063 39.339619 46.281905 29.311873 33.940063 40.882349 30.083238 59.395111 32.397333 31.625968 26.226413 58.623746 35.482793 37.796889 33.940063 20.826857
ENSRNOG00000000043 Idua protein_coding 30.218199 0.0641140 0.1209170 0.5302316 0.5959514 0.9992526 29.691417 23.7531337 31.670845 23.753134 33.650273 30.681131 27.711989 21.773706 23.753134 22.763420 57.403406 31.670845 39.588556 36.619414 27.711989 23.753134 28.701703 43.547412 55.423979 14.845709 25.732561 28.701703 30.681131 20.783992 31.670845 30.681131 24.742848 27.711989 42.557698 25.732561 47.506267 26.722275
ENSRNOG00000000044 Tmem175 protein_coding 24.418439 -0.0405021 0.1332847 -0.3038768 0.7612217 0.9992526 15.191208 18.9890099 15.191208 19.748570 17.469889 18.229450 19.748570 13.672087 12.152966 4.557362 32.661097 9.114725 17.469889 15.950768 22.027251 20.508131 21.267691 23.546372 31.901537 15.191208 18.229450 24.305933 26.584614 12.912527 29.622856 15.950768 23.546372 14.431648 23.546372 25.065493 21.267691 17.469889
ENSRNOG00000000047 Cd82 protein_coding 46.149014 0.0543979 0.1009968 0.5386105 0.5901556 0.9992526 31.675193 9.4351639 33.023074 22.240029 36.392775 37.066715 24.935790 27.631552 20.892149 21.566089 55.937043 16.174567 28.305492 28.305492 29.653372 24.261850 42.458238 50.545521 47.849760 28.979432 45.153999 33.697014 36.392775 45.153999 43.806118 32.349133 28.979432 21.566089 32.349133 45.153999 36.392775 33.023074
ENSRNOG00000000048 Gak protein_coding 94.890645 -0.0647427 0.0736958 -0.8785128 0.3796655 0.9992526 21.717392 18.4759901 33.062298 28.524336 28.524336 30.469177 18.475990 22.041532 23.338093 16.855289 50.565868 22.365672 35.655420 30.145036 31.117457 25.931214 36.303700 49.593447 53.807269 38.248541 27.227775 32.089878 41.165803 33.062298 39.545102 39.869242 23.013953 31.117457 46.352045 40.841662 32.089878 37.276120 
write.table(de2,file="de2.tsv",quote=FALSE,sep="\t")

ss3 %>% kbl() %>% kable_paper("hover", full_width = F)
UDI fraction Group RatID label trt
11 11 wholetissue S 11 wholetissue S 0
12 12 wholetissue S 12 wholetissue S 0
15 15 wholetissue S 15 wholetissue S 0
16 16 wholetissue S 16 wholetissue S 0
19 19 wholetissue S 19 wholetissue S 0
20 20 wholetissue S 20 wholetissue S 0
23 23 wholetissue S 23 wholetissue S 0
24 24 wholetissue S 24 wholetissue S 0
29 29 wholetissue S 29 wholetissue S 0
3 3 wholetissue S 3 wholetissue S 0
30 30 wholetissue S 30 wholetissue S 0
31 31 wholetissue S 31 wholetissue S 0
32 32 wholetissue S 32 wholetissue S 0
35 35 mito S 3 mito S 1
36 36 mito S 4 mito S 1
39 39 mito S 7 mito S 1
4 4 wholetissue S 4 wholetissue S 0
40 40 mito S 8 mito S 1
43 43 mito S 11 mito S 1
44 44 mito S 12 mito S 1
47 47 mito S 15 mito S 1
48 48 mito S 16 mito S 1
51 51 mito S 19 mito S 1
52 52 mito S 20 mito S 1
55 55 mito S 23 mito S 1
56 56 mito S 24 mito S 1
61 61 mito S 29 mito S 1
62 62 mito S 30 mito S 1
63 63 mito S 31 mito S 1
64 64 mito S 32 mito S 1
7 7 wholetissue S 7 wholetissue S 0
8 8 wholetissue S 8 wholetissue S 0 
de3 <- run_de(ss3,x3)
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 211 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
##  chr [1:21698] "ENSRNOG00000000500 Scube3 protein_coding" ...
##  chr [1:6568] "ENSRNOG00000000130 Dnajb5 protein_coding" ...

as.data.frame(de3[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)
baseMean log2FoldChange lfcSE stat pvalue padj 11 12 15 16 19 20 23 24 29 3 30 31 32 35 36 39 4 40 43 44 47 48 51 52 55 56 61 62 63 64 7 8
ENSRNOG00000000001 Arsj protein_coding 37.304995 3.4533016 0.1594917 21.651916 0.0000000 0.0000000 0.5326559 0.5681663 0.4261247 0.5326559 0.7102078 0.5681663 0.8522494 1.0298014 0.5681663 0.5681663 0.3551039 0.8522494 0.9232702 1.0653117 0.4616351 0.7812286 0.4616351 0.8522494 1.0298014 0.9587806 0.2840831 0.9232702 1.3138845 1.1008221 0.4971455 0.5681663 1.1363325 1.1718429 1.0653117 1.0653117 0.4971455 0.5326559
ENSRNOG00000000007 Gad1 protein_coding 53.747696 5.4076584 0.1734495 31.177140 0.0000000 0.0000000 0.3327388 0.1996433 0.1996433 0.0998216 0.6322037 0.1663694 0.2661910 0.0665478 0.1663694 0.3660127 0.2994649 0.2994649 0.1996433 1.3642291 1.1313119 0.9316687 0.1996433 0.5656560 1.4640507 1.1645858 1.3642291 1.0647642 1.7302418 1.5971463 0.8651209 1.2644075 1.4640507 1.7635157 1.5971463 1.3642291 0.1330955 0.0998216
ENSRNOG00000000008 Alx4 protein_coding 29.781602 2.3063794 0.2158722 10.684006 0.0000000 0.0000000 1.8348617 0.6382128 1.5556436 0.9174309 2.1539681 1.1966490 1.5157553 0.5584362 1.2365373 1.2764255 1.3562021 1.1168724 0.8775426 0.6781011 0.2393298 0.5185479 0.7179894 0.6382128 1.5955319 0.5983245 1.2365373 0.2393298 0.9972075 0.7179894 0.9972075 0.1994415 1.0769841 0.5983245 0.7179894 0.8376543 0.5584362 0.5584362
ENSRNOG00000000009 Tmco5b protein_coding 28.271135 9.7504099 0.5327924 18.300580 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.7203024 0.3961663 0.9363931 0.0000000 0.7923326 0.7203024 0.3601512 0.7563175 0.2521058 1.2965443 1.1884990 0.8643629 0.6122570 0.8283478 0.7563175 0.9003780 0.6482722 0.0000000 0.0000000
ENSRNOG00000000010 Cbln1 protein_coding 30.268375 6.9617446 0.3286441 21.183236 0.0000000 0.0000000 0.0324797 0.0324797 0.0000000 0.0324797 0.0324797 0.0649594 0.0974391 0.0974391 0.0974391 0.0324797 0.0000000 0.1299188 0.0000000 0.7795127 0.5196751 0.3897563 0.0000000 0.6820736 0.7795127 0.4222360 0.6820736 0.4222360 1.2017487 0.5196751 0.7795127 0.7795127 0.8769517 0.8444721 0.9743908 1.0393502 0.0000000 0.0649594
ENSRNOG00000000012 Tcf15 protein_coding 84.666617 -2.5799513 0.1776731 -14.520777 0.0000000 0.0000000 10.9843052 17.2283597 10.0362551 11.3439105 16.3783837 10.4285517 27.2646148 20.5955723 15.4303335 16.1168526 9.6766498 15.2995680 9.4478101 0.3269138 0.1961483 0.0653828 8.0747720 0.2288397 0.6211363 0.4576794 0.1961483 0.3596052 0.2615311 0.4903708 0.1961483 0.2288397 0.3269138 0.2942225 0.3269138 0.1307655 11.4746760 8.7939825
ENSRNOG00000000017 Steap1 protein_coding 14.467609 6.2530996 0.3885924 16.091667 0.0000000 0.0000000 0.0623921 0.0623921 0.0000000 0.0000000 0.0311960 0.0000000 0.0623921 0.1247841 0.0000000 0.0935881 0.0000000 0.0000000 0.0000000 0.4679405 0.3119604 0.1247841 0.0000000 0.2807643 0.5303326 0.3119604 0.2495683 0.2495683 0.4055485 0.2183723 0.3431564 0.4367445 0.3119604 0.1871762 0.3743524 0.4679405 0.0000000 0.0935881
ENSRNOG00000000021 AABR07061902.1 protein_coding 10.504602 8.3140537 0.5796121 14.344169 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.5239202 0.1746401 0.3492801 0.0000000 0.3492801 0.5937762 0.2444961 0.1047840 0.2095681 0.4191361 0.2095681 0.2095681 0.1397120 0.3842081 0.2095681 0.1746401 0.0000000 0.0000000 0.0000000
ENSRNOG00000000024 Hebp1 protein_coding 25.248485 0.6360841 0.1401613 4.538228 0.0000057 0.0000066 3.2086295 2.4270403 2.0156775 2.1390864 2.3447677 3.0440844 3.2909021 3.0440844 2.0979501 2.2213589 1.9745413 2.5093128 1.3163608 0.3702265 0.1645451 0.2879539 1.6043148 0.2879539 0.9461344 0.4524990 0.2879539 0.3702265 0.4113628 0.4524990 0.5759079 0.6170441 0.6170441 0.6170441 0.7404530 0.4524990 2.8384031 1.5220422
ENSRNOG00000000033 Tmcc2 protein_coding 223.839246 -1.7211422 0.1031026 -16.693497 0.0000000 0.0000000 32.5431565 49.9441245 37.1443740 40.9090065 42.3730302 48.1036375 34.9692530 56.0511950 35.2620577 38.7757147 46.8905892 32.9614490 20.7891372 1.5476822 1.5058530 1.0875605 46.3886382 1.4640237 2.1751210 1.4640237 1.1712190 1.0875605 2.5515842 1.5895115 1.4640237 1.3385360 1.9241455 1.5058530 2.3006087 1.7149992 44.9246145 42.1638840
ENSRNOG00000000034 Nuak2 protein_coding 58.253462 2.3100060 0.1188663 19.433646 0.0000000 0.0000000 2.7052276 2.4797919 2.0664933 1.8786303 2.0289207 1.5780494 3.8699783 3.6445427 1.7659124 1.5404768 1.2398960 2.4797919 1.4277590 1.1647508 0.8641699 0.9393151 1.6907672 0.6763069 2.0664933 1.1647508 1.0896056 0.8265973 2.1416385 1.8786303 1.5780494 1.3526138 1.5780494 2.0289207 1.3150412 1.4277590 2.4797919 1.2023234
ENSRNOG00000000035 Tmed11 protein_coding 5.675573 7.4293401 0.5638962 13.175013 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1850489 0.1850489 0.1480391 0.0000000 0.1480391 0.1480391 0.0370098 0.1110293 0.0740196 0.2590685 0.1850489 0.1850489 0.1480391 0.1850489 0.1480391 0.1480391 0.1850489 0.0000000 0.0000000
ENSRNOG00000000036 Klhdc8a protein_coding 45.871215 6.6164374 0.2618623 25.266854 0.0000000 0.0000000 0.0000000 0.1915987 0.0957993 0.0000000 0.0957993 0.0000000 0.0957993 0.0478997 0.1436990 0.2873980 0.1436990 0.1915987 0.0478997 1.1016924 1.3411907 1.2453914 0.0478997 0.6705954 3.6882745 1.6285887 1.3890904 0.8621940 2.6344818 1.4848897 1.5806891 1.7243881 1.7722878 2.0117861 2.0596858 1.3411907 0.0957993 0.5268964
ENSRNOG00000000040 RGD1304622 protein_coding 97.528259 1.7879746 0.0994641 17.976074 0.0000000 0.0000000 82.3432858 112.2862988 51.7197497 100.0368843 76.2185786 100.7174074 73.4964865 121.8136211 71.4549174 79.6211937 87.1069469 80.9822397 50.3587037 35.3871972 36.0677202 20.4156907 94.5927002 27.9014439 63.9691641 28.5819670 40.8313814 36.0677202 62.6081181 24.4988288 36.7482432 27.9014439 45.5950425 42.1924274 44.2339965 40.1508583 90.5095620 80.9822397
ENSRNOG00000000041 Slc26a1 protein_coding 46.021902 5.0993352 0.8887539 5.737623 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.8561578 1.7123156 0.0000000 0.0000000 0.0000000 0.0000000 2.5684734 0.8561578 0.0000000 0.0000000 40.2394160 18.8354713 18.8354713 0.0000000 20.5477869 34.2463115 17.9793135 32.5339959 23.1162602 54.7940983 23.1162602 29.1093647 19.6916291 34.2463115 36.8147848 29.1093647 25.6847336 24.8285758 0.0000000
ENSRNOG00000000042 Xpr1 protein_coding 136.048749 -0.6202520 0.0810839 -7.649512 0.0000000 0.0000000 293.2119743 298.8290619 232.5474279 262.8797011 312.3100722 272.9904588 274.1138763 329.7230439 218.5047088 236.4793892 271.3053325 221.3132526 237.0410980 21.3449330 20.7832242 14.6044278 261.7562836 12.9193015 29.2088557 19.0980979 15.7278454 13.4810103 33.7025258 21.3449330 21.9066418 43.2515748 25.8386031 27.5237294 24.7151856 15.1661366 251.0838171 222.4366701
ENSRNOG00000000043 Idua protein_coding 104.427192 -0.5102974 0.1039698 -4.908131 0.0000009 0.0000011 15.3736087 16.0319244 13.9020794 15.3348842 12.3918257 11.6560610 12.9339680 14.3667729 12.0820300 10.1845317 14.0569772 14.5603951 11.7722344 1.2391826 0.9293869 1.0842847 12.2369278 0.8519380 2.2460184 1.2391826 1.0842847 0.9293869 2.1685695 0.5808668 1.2004581 0.8132136 1.6651516 1.0068358 1.8587738 1.0455603 12.4305501 14.2118750
ENSRNOG00000000044 Tmem175 protein_coding 80.633148 -0.2929697 0.0950278 -3.082989 0.0020493 0.0022562 424.1173934 374.9647154 332.8338485 369.3472665 443.7784647 303.3422417 457.8220870 433.9479291 311.7684151 268.2331859 374.9647154 324.4076751 287.8942571 28.0872446 36.5134180 36.5134180 337.0469352 25.2785201 60.3875759 16.8523468 40.7265047 37.9177802 58.9832137 28.0872446 49.1526780 23.8741579 43.5352291 46.3439536 39.3221424 32.3003313 303.3422417 292.1073438
ENSRNOG00000000047 Cd82 protein_coding 90.873373 1.4571363 0.0810879 17.969842 0.0000000 0.0000000 91.1741644 112.3500994 70.5864499 67.6453478 74.1157724 84.1155194 117.6440831 103.5267931 69.9982295 92.9388257 71.7628907 76.4686540 58.8220416 28.8228004 19.4112737 21.7641554 69.4100090 24.1170370 48.8222945 14.1172900 25.8816983 21.1759350 41.7636495 25.2934779 31.7639024 39.4107678 28.2345799 39.4107678 31.7639024 28.8228004 82.9390786 62.3513641
ENSRNOG00000000048 Gak protein_coding 253.154440 0.3647060 0.0575672 6.335312 0.0000000 0.0000000 949.8535642 983.1817595 676.2846283 872.0877754 898.4725966 908.1933202 1073.4456216 1312.2976874 723.4995716 744.3296936 788.7672873 938.7441658 681.8393275 141.6448298 122.2033825 79.1544637 749.8843928 94.4298865 216.6332690 95.8185613 133.3127809 111.0939841 230.5200171 163.8636266 176.3616998 141.6448298 198.5804966 174.9730250 137.4788054 159.6976022 913.7480194 830.4275313 
write.table(de3,file="de3.tsv",quote=FALSE,sep="\t")

ss4 %>% kbl() %>% kable_paper("hover", full_width = F)
UDI fraction Group RatID label trt
1 1 wholetissue T 1 wholetissue T 0
10 10 wholetissue T 10 wholetissue T 0
13 13 wholetissue T 13 wholetissue T 0
14 14 wholetissue T 14 wholetissue T 0
17 17 wholetissue T 17 wholetissue T 0
18 18 wholetissue T 18 wholetissue T 0
2 2 wholetissue T 2 wholetissue T 0
21 21 wholetissue T 21 wholetissue T 0
22 22 wholetissue T 22 wholetissue T 0
25 25 wholetissue T 25 wholetissue T 0
26 26 wholetissue T 26 wholetissue T 0
27 27 wholetissue T 27 wholetissue T 0
28 28 wholetissue T 28 wholetissue T 0
33 33 mito T 1 mito T 1
34 34 mito T 2 mito T 1
37 37 mito T 5 mito T 1
38 38 mito T 6 mito T 1
41 41 mito T 9 mito T 1
42 42 mito T 10 mito T 1
45 45 mito T 13 mito T 1
46 46 mito T 14 mito T 1
49 49 mito T 17 mito T 1
5 5 wholetissue T 5 wholetissue T 0
50 50 mito T 18 mito T 1
53 53 mito T 21 mito T 1
54 54 mito T 22 mito T 1
57 57 mito T 25 mito T 1
58 58 mito T 26 mito T 1
59 59 mito T 27 mito T 1
6 6 wholetissue T 6 wholetissue T 0
60 60 mito T 28 mito T 1
9 9 wholetissue T 9 wholetissue T 0 
de4 <- run_de(ss4,x4)
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 72 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
##  chr [1:21715] "ENSRNOG00000000156 Megf6 protein_coding" ...
##  chr [1:6579] "ENSRNOG00000000383 Mypn protein_coding" ...

as.data.frame(de4[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)
baseMean log2FoldChange lfcSE stat pvalue padj 1 10 13 14 17 18 2 21 22 25 26 27 28 33 34 37 38 41 42 45 46 49 5 50 53 54 57 58 59 6 60 9
ENSRNOG00000000001 Arsj protein_coding 35.970770 3.1974984 0.1744971 18.324079 0.0000000 0.0000000 0.8052414 1.4315403 1.3868046 0.8052414 2.1920460 0.8947127 0.9841839 0.7157701 0.8052414 0.8947127 0.5368276 1.2078621 0.7157701 0.6262989 0.5368276 1.2525977 1.5657472 0.7157701 0.7605058 1.6104828 0.6262989 1.1183908 0.8499770 0.8499770 1.3868046 0.8947127 0.9394483 1.1631265 0.6262989 0.8052414 1.0736552 0.5815632
ENSRNOG00000000007 Gad1 protein_coding 55.872777 5.2656174 0.1554311 33.877510 0.0000000 0.0000000 0.2454429 0.1534018 0.3681643 0.0306804 0.2761232 0.2454429 0.0920411 0.6442876 0.1534018 0.3681643 0.1534018 0.3374840 0.3068036 1.0431323 0.6442876 1.1351733 1.2885751 0.7670090 0.5522465 1.2578948 1.1658537 1.2272144 0.2454429 1.5033377 1.0124519 1.2885751 1.7794609 1.1658537 1.0431323 0.5829269 1.5033377 0.1840822
ENSRNOG00000000008 Alx4 protein_coding 26.988602 1.9097179 0.2097035 9.106752 0.0000000 0.0000000 0.2000711 1.5338784 0.8336295 1.4671880 2.4341983 1.3338073 0.7002488 2.1340916 0.3667970 1.1670814 1.4671880 2.0674013 0.9670103 0.4668325 0.3334518 0.5668681 0.8336295 0.2000711 0.2667615 0.5668681 0.3001066 0.6335585 0.8336295 0.5668681 0.4334874 0.4668325 0.6669036 0.4001422 0.4001422 0.9336651 0.7335940 0.8336295
ENSRNOG00000000009 Tmco5b protein_coding 26.299815 8.9045951 0.5349989 16.644137 0.0000000 0.0000000 0.0360027 0.0360027 0.0000000 0.0360027 0.0000000 0.0000000 0.0000000 0.0720055 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.6480494 0.2520192 0.8280631 1.0440796 0.1800137 0.2880219 0.7920604 0.4680357 0.7920604 0.0000000 0.9720741 0.7920604 1.1880905 1.0080768 0.5400412 0.5400412 0.0000000 0.3600274 0.0000000
ENSRNOG00000000010 Cbln1 protein_coding 32.075967 7.7195551 0.3915108 19.717349 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0328934 0.0328934 0.0000000 0.1315738 0.0328934 0.1973607 0.0000000 0.0000000 0.0000000 0.0657869 0.7565493 0.3618279 0.9868034 0.5920820 0.5591886 0.4605082 0.6578689 0.8223361 0.5920820 0.0000000 0.8223361 0.7894427 0.7565493 0.5591886 0.8881230 0.8881230 0.0000000 0.8223361 0.0000000
ENSRNOG00000000012 Tcf15 protein_coding 110.022720 -2.6384953 0.1656088 -15.932099 0.0000000 0.0000000 8.2486227 20.1171447 20.8292560 13.2334020 10.9190402 14.3015689 9.5541602 26.8822022 17.7434403 18.3665377 19.2270055 25.7250213 22.4908491 0.5637548 0.3263844 0.1483565 0.5340835 0.1186852 0.1780278 0.3857270 0.2967130 0.1186852 24.6865256 0.4450696 0.2967130 0.4450696 0.3560557 0.2373704 0.1483565 19.1676629 0.2967130 12.9960315
ENSRNOG00000000017 Steap1 protein_coding 12.370744 7.7704079 0.5621036 13.823801 0.0000000 0.0000000 0.0000000 0.0000000 0.0513400 0.0000000 0.1540201 0.0000000 0.0513400 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.5647404 0.2053601 1.0268007 0.3080402 0.3593803 0.3593803 0.6160804 0.3080402 0.6160804 0.0000000 0.3080402 0.3593803 0.2053601 0.5647404 0.1540201 0.3080402 0.0000000 0.4620603 0.0000000
ENSRNOG00000000021 AABR07061902.1 protein_coding 11.716808 8.4363830 0.5602052 15.059452 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0253799 0.0000000 0.0000000 0.0761396 0.1015194 0.3806979 0.2030389 0.2030389 0.1015194 0.3045583 0.2030389 0.2284187 0.0000000 0.1776590 0.0761396 0.2284187 0.3299382 0.3806979 0.2030389 0.0000000 0.0507597 0.0000000
ENSRNOG00000000024 Hebp1 protein_coding 22.989227 0.3159324 0.1566626 2.016642 0.0437329 0.0464945 1.4623450 3.3633934 3.1074830 1.9010485 2.7784554 2.2300761 1.3161105 2.7784554 2.7418968 2.0838416 2.3763106 2.8515727 2.6687796 0.1827931 0.2193517 0.2193517 0.2559104 0.3290276 0.1096759 0.4752621 0.3290276 0.1827931 2.0107243 0.4021449 0.2193517 0.5849380 0.4021449 0.4021449 0.3655862 2.4859864 0.4021449 1.4257863
ENSRNOG00000000033 Tmcc2 protein_coding 224.124493 -1.7154163 0.0849979 -20.181865 0.0000000 0.0000000 30.8389182 29.3246370 33.2850647 28.8295836 30.8097974 34.7411042 18.5790650 35.7020903 31.5960587 39.4877932 35.5273656 38.9927398 38.4103239 0.9027445 0.6697782 1.3977980 1.1939525 0.5532950 0.9318653 1.1065901 1.1357109 0.9609861 27.3735440 1.3395564 1.1357109 1.1648317 1.2521940 1.2230732 0.9027445 27.2861816 0.6406574 18.6373065
ENSRNOG00000000034 Nuak2 protein_coding 61.810761 2.2596789 0.1269558 17.798942 0.0000000 0.0000000 2.1659042 2.9470500 3.0180633 1.8818512 2.2014108 1.3492518 1.3847584 2.1303976 2.6274904 3.0890765 2.4854639 4.0477554 1.4557717 1.2072253 0.8521590 0.9586789 1.9173578 0.7811458 0.6036127 1.2427319 1.1717187 1.1717187 3.2666096 1.5977982 0.8876657 1.6333048 1.4202651 1.6688114 1.4202651 1.9173578 1.3492518 1.9173578
ENSRNOG00000000035 Tmed11 protein_coding 6.960863 7.8382465 0.5684970 13.787667 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1637000 0.1964400 0.0982200 0.0654800 0.0654800 0.1309600 0.1964400 0.0982200 0.1964400 0.0000000 0.3274000 0.1309600 0.1637000 0.2291800 0.0654800 0.1309600 0.0000000 0.1637000 0.0000000
ENSRNOG00000000036 Klhdc8a protein_coding 45.197897 6.8182015 0.2539232 26.851427 0.0000000 0.0000000 0.0323555 0.2264887 0.0323555 0.1941331 0.2588442 0.0000000 0.0000000 0.0647110 0.0970666 0.1617776 0.0647110 0.0647110 0.0323555 1.0030212 0.4529773 0.7441770 1.4236430 0.7441770 0.4206218 0.9706657 0.8412436 1.2618654 0.0000000 1.1324433 1.4236430 1.1000878 1.3265765 1.1971544 1.0353768 0.0000000 0.9059547 0.0647110
ENSRNOG00000000040 RGD1304622 protein_coding 92.744370 1.7260883 0.0903456 19.105390 0.0000000 0.0000000 21.0083343 30.8731173 28.4982622 22.2871025 23.0178271 27.2194940 15.1625369 27.2194940 21.0083343 22.8351460 24.6619576 27.2194940 24.6619576 10.9608701 4.0189857 8.2206525 8.7686960 6.7592032 6.5765220 10.4128266 7.1245655 9.8647831 30.1423927 8.5860149 4.7497104 11.1435512 13.5184064 9.8647831 11.3262324 26.3060881 9.1340584 18.4507979
ENSRNOG00000000041 Slc26a1 protein_coding 43.746306 7.1829272 0.4404495 16.308173 0.0000000 0.0000000 1.1013186 0.1001199 0.0000000 0.1001199 0.0000000 0.1001199 0.8009590 0.0000000 0.3003596 0.0000000 0.0000000 0.1001199 0.0000000 3.2038360 0.7008391 3.1037161 3.5041956 2.3027571 2.2026372 3.9046751 2.7032366 2.3027571 0.0000000 4.3051546 3.5041956 2.5029968 3.7044353 4.3051546 2.5029968 0.1001199 3.0035962 0.0000000
ENSRNOG00000000042 Xpr1 protein_coding 136.192281 -0.4154707 0.0772851 -5.375817 0.0000001 0.0000001 342.0182326 378.9532469 372.3049443 333.1538292 372.3049443 386.3402498 241.5549937 489.0195896 373.0436446 415.1495610 378.2145466 375.9984458 405.5464572 16.9901066 18.4675072 31.7641123 31.7641123 20.6836080 16.2514063 25.8545100 24.3771095 32.5028126 302.8671174 37.6737146 32.5028126 39.1511152 31.0254120 30.2867117 25.1158097 310.9928206 56.1412218 259.2838005
ENSRNOG00000000043 Idua protein_coding 106.370254 -0.4333749 0.0795492 -5.447883 0.0000001 0.0000001 325.7140296 412.0610375 473.2690431 367.2480334 377.0850343 440.4790401 311.5050283 463.4320422 417.5260380 489.6640446 400.0380364 509.3380463 488.5710445 32.7900030 26.2320024 37.1620034 33.8830031 26.2320024 25.1390023 43.7200040 40.4410037 31.6970029 436.1070397 48.0920044 28.4180026 31.6970029 34.9760032 33.8830031 27.3250025 380.3640346 30.6040028 301.6680275
ENSRNOG00000000044 Tmem175 protein_coding 81.215109 -0.3851394 0.0923749 -4.169308 0.0000306 0.0000351 390.3986021 552.0614594 490.1480247 460.9111250 462.6309426 481.5489366 297.5284501 708.5648638 459.1913074 660.4099701 433.3940429 605.3758059 581.2983591 34.3963526 42.9954408 39.5558055 41.2756231 27.5170821 10.3189058 39.5558055 36.1161702 48.1548937 460.9111250 53.3143465 41.2756231 55.0341642 67.0728876 36.1161702 53.3143465 447.1525839 32.6765350 421.3553195
ENSRNOG00000000047 Cd82 protein_coding 93.529245 1.4943004 0.0893353 16.726871 0.0000000 0.0000000 235.6725245 307.3169720 280.9216493 233.7871443 380.8467997 284.6924096 201.7356810 431.7520650 279.0362691 280.9216493 258.2970869 341.2538155 309.2023522 88.6128692 26.3953227 101.8105306 103.6959108 58.4467861 60.3321663 79.1859682 79.1859682 118.7789524 352.5660967 141.4035147 126.3204732 94.2690098 122.5497128 90.4982494 81.0713484 246.9848057 60.3321663 218.7041028
ENSRNOG00000000048 Gak protein_coding 243.747442 0.3480978 0.0599274 5.808657 0.0000000 0.0000000 523.6878738 730.4879968 782.9596699 560.7267018 646.1217775 849.8353315 386.8499814 846.7487625 766.4979685 842.6333372 732.5457095 771.6422502 805.5945092 68.9333744 58.6448110 90.5393574 96.7124954 74.0776560 53.5005293 113.1741967 95.6836390 115.2319094 686.2471745 157.4150190 86.4239320 101.8567770 125.5204727 126.5493290 73.0487997 633.7755015 98.7702080 515.4570231 
write.table(de4,file="de4.tsv",quote=FALSE,sep="\t")

Second pass DE analysis

Here we are going to exclude the mito samples with < 40% mito reads. Contrast 1 is not repeated.

  1. Mito fraction S versus T: 0 DEGs

  2. Whole versus mito fraction in S: 20 DEGs

  3. Whole versus mito fraction in T: 0 DEGs

  4. Mito fraction S versus T taking into consideration the mitofrac: 1 DEG

ss0 <- subset(ss, fraction=="mito" & mtfrac>0.4)
ss <- subset(ss, fraction!="mito")
ss <- rbind(ss,ss0)

ss2 <- subset(ss,fraction=="mito")
ss2$trt <- grepl("T",ss2$Group)*1
x2 <- x[,which(colnames(x) %in% rownames(ss2))]

ss3 <- subset(ss,Group=="S")
ss3$trt <- grepl("mito",ss3$fraction)*1
x3 <- x[,which(colnames(x) %in% rownames(ss3))]

ss4 <- subset(ss,Group=="T")
ss4$trt <- grepl("mito",ss4$fraction)*1
x4 <- x[,which(colnames(x) %in% rownames(ss4))]

ss2 %>% kbl() %>% kable_paper("hover", full_width = F)
UDI fraction Group RatID label nreads mtfrac trt
33 33 mito T 1 mito T 5474018 0.8850587 1
34 34 mito T 2 mito T 9988027 0.9542655 1
35 35 mito S 3 mito S 1469458 0.4723449 0
36 36 mito S 4 mito S 1168009 0.5294668 0
37 37 mito T 5 mito T 1353729 0.4835828 1
39 39 mito S 7 mito S 1780282 0.6982787 0
43 43 mito S 11 mito S 1700043 0.4096285 0
46 46 mito T 14 mito T 1086442 0.4345745 1
47 47 mito S 15 mito S 8344505 0.9222841 0
48 48 mito S 16 mito S 1296403 0.5776151 0
52 52 mito S 20 mito S 3085085 0.7437176 0
53 53 mito T 21 mito T 1628001 0.5591495 1
56 56 mito S 24 mito S 2133695 0.6219164 0
57 57 mito T 25 mito T 1841975 0.5667390 1
58 58 mito T 26 mito T 3802593 0.7906639 1
59 59 mito T 27 mito T 7532868 0.9057456 1
60 60 mito T 28 mito T 6001091 0.8726930 1
61 61 mito S 29 mito S 1785865 0.5147013 0
62 62 mito S 30 mito S 1474773 0.4122207 0
63 63 mito S 31 mito S 1503193 0.4336855 0
64 64 mito S 32 mito S 4091353 0.8123151 0 
de2 <- run_de(ss2,x2)
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
##  chr(0) 
##  chr(0)

as.data.frame(de2[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)
baseMean log2FoldChange lfcSE stat pvalue padj 33 34 35 36 37 39 43 46 47 48 52 53 56 57 58 59 60 61 62 63 64
ENSRNOG00000000001 Arsj protein_coding 21.954498 -0.1611982 0.1860877 -0.8662486 0.3863539 0.9770360 2.5575363 6.7405052 9.7242053 1.7257703 6.8437018 16.2514063 22.3695872 3.6816983 5.322005 22.2601024 3.7150196 16.8297615 10.8491273 14.290983 23.931328 6.561388 13.438866 3.203836 19.411274 18.4275071 4.9990910
ENSRNOG00000000007 Gad1 protein_coding 38.259273 0.0157898 0.1278073 0.1235441 0.9016762 0.9989235 3.4040757 12.3526287 25.1842597 5.6656365 6.7592032 15.7278454 14.2621678 5.0445594 10.021135 23.6384092 37.0255237 8.6782887 25.2795217 49.657152 4.553895 18.458448 33.225452 29.943013 48.783092 22.4961862 22.9580623
ENSRNOG00000000008 Alx4 protein_coding 15.586120 -0.1843950 0.2540039 -0.7259532 0.4678675 0.9821127 9.5273223 9.2043570 7.9673993 3.3597164 1.7020378 7.6468654 24.5700095 1.4997273 5.663116 3.3702526 5.8345232 1.7257703 1.2220896 14.774006 9.256381 3.155741 14.635513 23.116260 1.797590 9.7721196 14.2394796
ENSRNOG00000000009 Tmco5b protein_coding 17.958991 -0.1490540 0.2098614 -0.7102501 0.4775491 0.9828999 15.4108402 0.8388754 10.8579107 7.4587750 15.6520295 23.9313281 9.3734109 7.2793856 2.102517 4.1175429 20.2702578 3.6660001 3.1055799 15.727845 4.862103 1.991273 2.444179 16.990107 16.198667 6.5744612 11.9745103
ENSRNOG00000000010 Cbln1 protein_coding 21.660571 0.2118931 0.1682983 1.2590325 0.2080186 0.9691549 16.9901066 8.4850158 6.3114827 10.6440091 25.6847336 1.4380721 13.0294928 16.9517614 14.290983 11.9656641 7.4987287 13.4388658 2.4028770 9.999747 16.584756 4.499182 4.567029 15.166137 8.427645 3.9825469 7.8213735
ENSRNOG00000000012 Tcf15 protein_coding 9.262726 0.3178989 0.2847919 1.1162498 0.2643152 0.9712115 10.6724665 3.5655420 1.3275156 1.4665075 3.6935014 1.5427302 4.9965905 6.6525057 5.136947 1.3182328 8.1434330 6.7807046 4.7636612 11.045228 3.749364 2.799764 1.001199 5.882204 5.528252 1.6663637 0.7307247
ENSRNOG00000000017 Steap1 protein_coding 9.219725 -0.1566629 0.2798756 -0.5597590 0.5756438 0.9947511 6.4704246 2.4570009 2.4995455 1.8268117 11.2341753 1.2965607 2.2567766 1.4665075 5.909602 6.1709206 1.8408491 4.6567540 11.9862090 1.318233 1.628687 4.068423 6.124707 9.204357 2.812023 6.7194329 1.5017981
ENSRNOG00000000021 AABR07061902.1 protein_coding 7.215800 0.2243939 0.3701817 0.6061723 0.5444004 0.9931259 2.7613071 1.8746822 8.3992911 0.5005994 8.8233062 6.1425024 2.8328182 1.4614493 1.685126 1.9448411 0.7965094 0.7332538 2.9548011 10.027746 3.944677 5.322005 1.712316 1.318233 3.257373 3.3903523 0.0000000
ENSRNOG00000000024 Hebp1 protein_coding 9.716458 -0.3710322 0.2496551 -1.4861792 0.1372317 0.9527937 0.5991967 3.2573732 6.1026341 2.7220921 5.5226142 3.2806938 12.8789130 0.9010789 4.117543 5.5282521 1.8330000 1.0960870 8.4256314 3.565542 1.460267 2.444179 8.125703 11.570476 3.944677 11.9745103 9.4177357
ENSRNOG00000000033 Tmcc2 protein_coding 35.369665 0.0426072 0.1377696 0.3092642 0.7571206 0.9989235 23.9123174 6.0485043 24.6142711 30.8216803 5.7522885 14.1152839 35.2596637 26.5403979 25.772199 12.1854342 21.2782041 3.9046751 18.8230533 26.412760 6.998727 5.663116 12.357593 14.910448 4.779056 13.4429857 30.2867117
ENSRNOG00000000034 Nuak2 protein_coding 34.038444 0.0460973 0.1384433 0.3329689 0.7391577 0.9989235 11.0207660 3.1860375 7.5769556 16.9901066 20.8268571 6.5744612 36.5887813 28.2532070 3.475341 11.9437017 33.9035228 17.0130756 33.1356851 18.746822 26.317779 4.004795 22.352376 25.798510 8.998364 6.3938409 21.3449330
ENSRNOG00000000035 Tmed11 protein_coding 4.232297 0.2003516 0.3653643 0.5483613 0.5834438 0.9959132 3.0712512 0.9998182 0.9134058 2.8085438 0.9724205 0.5310063 0.9776717 2.2161009 2.314095 0.5259569 3.3262528 3.4246311 0.4793574 3.800269 1.356141 2.722092 4.602178 2.343353 2.239811 0.4004795 2.9411021
ENSRNOG00000000036 Klhdc8a protein_coding 31.161089 -0.0394331 0.1598930 -0.2466218 0.8052009 0.9989235 14.5287869 7.8393384 2.3027571 16.4701716 14.1277554 4.3325455 14.0664499 14.6044278 9.400065 2.3895281 7.5769556 32.5028126 27.7691428 10.782116 24.614271 27.397049 3.355502 20.087135 28.478959 29.2624900 25.7721995
ENSRNOG00000000040 RGD1304622 protein_coding 50.406120 -0.1011744 0.1538829 -0.6574764 0.5108747 0.9847052 32.5737320 14.9175500 35.3871972 48.7830920 21.0901746 16.7985822 9.4112681 22.9405962 36.855014 8.8317274 6.5765220 14.6044278 13.2897473 9.823616 13.198568 45.799418 38.568254 17.619556 41.245535 55.6502561 7.0705213
ENSRNOG00000000041 Slc26a1 protein_coding 29.927087 -0.0276829 0.1548625 -0.1787579 0.8581278 0.9989235 8.4153103 4.6567540 40.2394160 2.6364656 16.8297615 14.9175500 27.2209209 24.8517638 17.809481 15.1187240 2.7032366 20.5877145 14.1277554 6.165546 7.855290 14.042719 9.724205 5.310062 10.509971 25.1158097 23.1409523
ENSRNOG00000000042 Xpr1 protein_coding 38.502690 0.1192935 0.1623686 0.7347079 0.4625174 0.9820810 3.0532859 6.1104481 28.0706109 28.5405079 11.3080732 17.2965148 44.5202049 3.9546983 15.201075 16.2736909 25.8598749 40.4991707 36.0876320 23.518015 4.104915 19.999494 46.683018 7.665273 8.951377 24.7151856 8.7517848
ENSRNOG00000000043 Idua protein_coding 29.604229 0.0817579 0.1676953 0.4875384 0.6258769 0.9960496 4.9990910 4.3843480 17.9746804 7.7793643 4.5135531 6.8437018 42.8446166 28.5405079 7.363397 15.9660137 12.8423668 3.1158229 11.4008062 21.698255 21.096214 23.010892 13.122775 24.077968 2.603117 28.2345799 16.5847564
ENSRNOG00000000044 Tmem175 protein_coding 25.636539 0.0135423 0.1760225 0.0769350 0.9386753 0.9989235 11.1990548 2.5029968 11.7644083 15.9705062 3.8326364 4.7497104 24.1534768 6.8069437 3.849795 6.5992839 14.7740057 18.5127618 4.4706336 25.944772 17.979313 3.715020 10.315015 21.020184 22.457260 25.7721995 10.7794226
ENSRNOG00000000047 Cd82 protein_coding 46.089148 -0.0218677 0.1348713 -0.1621376 0.8711975 0.9989235 31.8693114 9.5273223 45.1013492 15.4661280 30.2374480 3.7044353 48.8222945 25.7985100 7.332000 6.5765220 24.1534768 21.7173919 8.8943547 15.887165 35.457614 33.168698 8.415310 31.932027 57.362572 6.4713245 26.6018811
ENSRNOG00000000048 Gak protein_coding 94.626970 -0.1295233 0.0951233 -1.3616358 0.1733129 0.9644552 44.5717882 48.8009938 12.2236130 47.7748069 59.6702001 38.7898123 143.5879688 43.5863607 53.755463 8.0095899 69.4100090 51.5970199 16.9969094 22.287102 69.090178 23.013953 12.744150 34.951763 93.076236 76.3651426 30.2425214 
write.table(de2,file="de2b.tsv",quote=FALSE,sep="\t")

ss3 %>% kbl() %>% kable_paper("hover", full_width = F)
UDI fraction Group RatID label nreads mtfrac trt
11 11 wholetissue S 11 wholetissue S 28160771 0.3290185 0
12 12 wholetissue S 12 wholetissue S 30053603 0.2640899 0
15 15 wholetissue S 15 wholetissue S 25070009 0.3071025 0
16 16 wholetissue S 16 wholetissue S 27766116 0.3024493 0
19 19 wholetissue S 19 wholetissue S 30788467 0.3233271 0
20 20 wholetissue S 20 wholetissue S 30589099 0.3209945 0
23 23 wholetissue S 23 wholetissue S 32055355 0.3547227 0
24 24 wholetissue S 24 wholetissue S 28630316 0.2413714 0
29 29 wholetissue S 29 wholetissue S 24309444 0.3591212 0
3 3 wholetissue S 3 wholetissue S 23906716 0.3067514 0
30 30 wholetissue S 30 wholetissue S 26615136 0.2724498 0
31 31 wholetissue S 31 wholetissue S 27019886 0.3536612 0
32 32 wholetissue S 32 wholetissue S 20876971 0.3177831 0
4 4 wholetissue S 4 wholetissue S 25823475 0.2520651 0
7 7 wholetissue S 7 wholetissue S 25564135 0.2969021 0
8 8 wholetissue S 8 wholetissue S 22385138 0.2206553 0
35 35 mito S 3 mito S 1469458 0.4723449 1
36 36 mito S 4 mito S 1168009 0.5294668 1
39 39 mito S 7 mito S 1780282 0.6982787 1
43 43 mito S 11 mito S 1700043 0.4096285 1
47 47 mito S 15 mito S 8344505 0.9222841 1
48 48 mito S 16 mito S 1296403 0.5776151 1
52 52 mito S 20 mito S 3085085 0.7437176 1
56 56 mito S 24 mito S 2133695 0.6219164 1
61 61 mito S 29 mito S 1785865 0.5147013 1
62 62 mito S 30 mito S 1474773 0.4122207 1
63 63 mito S 31 mito S 1503193 0.4336855 1
64 64 mito S 32 mito S 4091353 0.8123151 1 
de3 <- run_de(ss3,x3)
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 547 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
##  chr [1:8730] "ENSRNOG00000065610 na lncRNA" "ENSRNOG00000066930 na lncRNA" ...
##  chr [1:4758] "ENSRNOG00000028837 AC128207.1 protein_coding" ...

as.data.frame(de3[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)
baseMean log2FoldChange lfcSE stat pvalue padj 11 12 15 16 19 20 23 24 29 3 30 31 32 35 36 39 4 43 47 48 52 56 61 62 63 64 7 8
ENSRNOG00000000001 Arsj protein_coding 37.382998 1.6058722 0.5519473 2.9094664 0.0036205 0.0150348 0.5326559 0.7663947 7.9830068 0.6170441 6.4828036 0.5196751 0.9293869 1.0298014 0.7663947 10.6440091 0.4113628 7.7793643 0.8444721 1.1617337 0.4616351 1.0537927 8.6482574 1.1929520 2.5931214 0.8444721 1.2004581 0.5681663 1.532789 21.9532688 1.2340883 9.7242053 0.4547157 0.5808668
ENSRNOG00000000007 Gad1 protein_coding 56.810911 1.9271602 0.7899382 2.4396343 0.0147021 0.0364680 0.3327388 4.0831381 1.4665075 0.1254877 8.9047404 0.1634569 4.7057633 0.0665478 3.4026151 2.6885971 0.3764632 4.2180349 0.1961483 24.1170370 1.1313119 19.0546446 1.4665075 1.8404870 19.2154924 1.0461243 28.2345799 1.2644075 29.943013 12.9541499 2.0078040 19.2154924 0.1307655 1.7646612
ENSRNOG00000000008 Alx4 protein_coding 30.190271 1.3646831 0.4286220 3.1838850 0.0014531 0.0090983 1.8348617 13.6985246 1.5255748 0.8641699 30.2374480 0.9358811 4.5538950 0.5584362 26.5408914 1.2517537 1.2774686 15.6786767 0.6863128 2.0372688 0.2393298 11.1300512 0.7041114 1.5029042 17.3585349 0.1871762 2.1571082 0.1994415 23.116260 0.5867595 0.6763069 11.7590075 0.4367445 1.6777508
ENSRNOG00000000009 Tmco5b protein_coding 27.186272 1.7606833 1.7173770 1.0252165 0.3052610 0.3429545 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 15.4273015 0.3961663 14.6044278 0.0000000 0.7401956 14.2394796 0.2444961 25.4550475 0.6122570 12.919301 0.9381224 0.9252445 12.2052682 0.0000000 0.0000000
ENSRNOG00000000010 Cbln1 protein_coding 29.492568 2.1595541 1.0377306 2.0810353 0.0374307 0.0736637 0.0324797 0.0387245 0.0000000 0.0478997 0.6652506 0.0822726 0.9724205 0.0974391 0.1161734 0.0355104 0.0000000 2.6610023 0.0000000 7.7793643 0.5196751 0.4646935 0.0000000 1.1495921 13.9702620 0.5347716 5.1862428 0.7795127 1.045560 0.9232702 1.4369901 21.2880182 0.0000000 0.6482804
ENSRNOG00000000012 Tcf15 protein_coding 107.992711 -1.3915410 0.4060606 -3.4269297 0.0006104 0.0058963 10.9843052 309.9921590 10.2150814 236.1414889 122.4533791 13.3435307 390.8712351 20.5955723 277.6400362 16.4040232 201.4348147 114.3875877 12.0886532 4.6867055 0.1961483 1.1764408 8.2186485 12.9299374 1.4665075 0.4601217 7.0300582 0.2288397 5.882204 0.2994649 6.8052302 0.9776717 14.6820667 126.0723768
ENSRNOG00000000017 Steap1 protein_coding 14.129584 1.9332698 1.1110711 1.7400054 0.0818581 0.1351054 0.0623921 0.2396787 0.0000000 0.0000000 0.0391173 0.0000000 1.1199055 0.1247841 0.0000000 0.1196649 0.0000000 0.0000000 0.0000000 8.3992911 0.3119604 0.4793574 0.0000000 14.5546824 0.3129384 0.3005808 3.9196692 0.4367445 1.198393 0.2393298 10.2738934 0.5867595 0.0000000 1.6798582
ENSRNOG00000000021 AABR07061902.1 protein_coding 9.742919 1.0559034 1.7466523 0.6045298 0.5454915 0.5674987 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 10.1710568 0.1746401 7.7136508 0.0000000 9.5490490 0.1340175 0.2220587 4.0684227 0.1397120 8.485016 0.2160907 2.8085438 0.0000000 0.0000000 0.0000000
ENSRNOG00000000024 Hebp1 protein_coding 28.538572 0.4621016 0.1739104 2.6571248 0.0078810 0.0238609 3.2086295 19.1242705 1.5915050 2.0136717 2.0240923 3.5445755 53.2200456 3.0440844 16.5311491 1.7539035 1.8587738 2.1661339 1.5327894 5.9872551 0.1645451 2.2689812 1.2667081 0.8906625 0.2485727 0.4310970 7.3177563 0.6170441 4.862103 0.4871954 0.6970402 0.3906143 3.3050772 24.6142711
ENSRNOG00000000033 Tmcc2 protein_coding 277.956987 -0.6604479 0.3106153 -2.1262567 0.0334819 0.0676945 32.5431565 559.5926316 29.0299495 575.2795665 33.7064411 782.6014762 204.3333831 56.0511950 395.0892700 30.3049135 659.3950859 26.2198180 338.2199423 9.0434631 1.5058530 12.1854342 36.2547455 30.5874616 0.9316687 17.6935986 9.2878811 1.3385360 21.558845 1.1768898 32.3521229 1.3642291 730.8817265 246.3732658
ENSRNOG00000000034 Nuak2 protein_coding 61.087509 1.3861424 0.3695315 3.7510801 0.0001761 0.0031598 2.7052276 36.9568808 1.7157820 5.9919672 2.1539681 35.9586270 4.0290821 3.6445427 26.3177788 1.2790375 3.9546983 2.6326277 32.5339959 1.2126364 0.8641699 13.9988185 1.4038216 6.5911639 1.1567607 18.8354713 1.9558651 1.3526138 23.518015 1.6845859 4.1943770 1.5157553 56.5064139 1.2517537
ENSRNOG00000000035 Tmed11 protein_coding 5.671059 1.5646015 1.5262539 1.0251253 0.3053041 0.3429772 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.2233625 0.1850489 2.7122818 0.0000000 3.0854603 0.1080454 1.1234175 0.2233625 0.1480391 3.390352 0.1397120 3.0854603 0.1800756 0.0000000 0.0000000
ENSRNOG00000000036 Klhdc8a protein_coding 45.974588 2.1359588 0.9808878 2.1775772 0.0294375 0.0613353 0.0000000 2.6610023 0.0822726 0.0000000 0.0649594 0.0000000 0.0710208 0.0478997 1.9957517 0.2468177 0.9724205 0.1299188 0.0387245 0.8167390 1.3411907 17.2965148 0.0411363 24.9587937 0.9419111 0.6970402 1.1008221 1.7243881 24.614271 1.7277236 13.9380276 0.9094314 0.0774489 0.3906143
ENSRNOG00000000040 RGD1304622 protein_coding 103.940290 1.0913142 0.2965219 3.6803833 0.0002329 0.0036382 82.3432858 40.3289572 3.1790230 68.8945702 3.6614351 87.0566215 3.5935791 121.8136211 25.6638819 4.8940222 59.9898298 3.8902748 43.5283108 1.7302418 36.0677202 7.3325377 5.8142657 44.0550313 1.9614831 31.1756820 1.1978597 27.9014439 16.376001 2.5934135 30.4635855 1.9287917 78.2333153 3.9595918
ENSRNOG00000000041 Slc26a1 protein_coding 44.439600 1.9099061 1.3612930 1.4030088 0.1606142 0.2193057 0.0000000 0.0000000 0.0000000 0.5599527 0.0623921 0.0000000 0.0000000 0.0000000 0.0000000 0.1127178 0.5599527 0.0000000 0.0000000 1.8747500 18.8354713 0.8605807 0.0000000 22.3981096 1.1854494 3.2356623 1.0769841 19.6916291 1.564692 1.6156220 19.0383932 0.9358811 3.4753410 0.0000000
ENSRNOG00000000042 Xpr1 protein_coding 163.194639 -0.2981857 0.1402995 -2.1253510 0.0335573 0.0677922 293.2119743 23.7657682 15.3220484 317.3369732 19.4199743 374.8834275 17.5753786 329.7230439 17.3776012 15.5811168 327.5080300 13.7616364 325.5160625 1.3685746 20.7832242 1.1614849 17.2465569 35.2596637 0.9779843 18.5127618 1.3685746 43.2515748 2.054935 1.8134792 29.8351000 0.9430563 344.8001894 14.2619875
ENSRNOG00000000043 Idua protein_coding 122.810154 -0.2686653 0.1535189 -1.7500474 0.0801101 0.1328574 15.3736087 14.7013020 17.1959812 263.4392257 13.1636083 97.5661935 10.8482179 14.3667729 11.0792421 12.5976129 241.4859569 15.4672398 98.5386140 1.0393502 0.9293869 0.9942910 15.1362954 38.5845331 1.1518157 7.7793643 0.4871954 0.8132136 1.526947 1.2453914 31.9320274 1.1106795 104.0489970 11.9200479
ENSRNOG00000000044 Tmem175 protein_coding 95.444315 -0.1834256 0.1051493 -1.7444303 0.0810841 0.1341243 177.6425655 8.8841261 161.2839564 64.2819136 13.2180430 101.2328379 10.6573914 181.7601084 7.3868015 129.9798973 65.2595853 9.6625567 96.0774619 0.6538277 15.2937308 0.8651209 163.3255255 10.5099707 1.2130482 12.6541047 0.6538277 9.9997471 1.031490 22.4572598 6.8437018 0.9620727 101.2328379 6.7998080
ENSRNOG00000000047 Cd82 protein_coding 99.321535 0.8880817 0.2440332 3.6391832 0.0002735 0.0039336 18.5750982 7.6186650 102.7389344 4.4984898 4.7341483 80.0732418 6.2392071 21.0917244 4.7467075 135.2729303 4.7723109 4.8844387 55.9952740 1.5286058 3.9546983 1.4758670 101.0266188 3.2467361 1.6531946 20.1582986 1.3414295 8.0292360 1.914638 57.3625717 2.1123343 1.8410577 78.9533363 3.3067798
ENSRNOG00000000048 Gak protein_coding 288.809504 0.2798440 0.0760348 3.6804724 0.0002328 0.0036382 527.6137127 25.4987050 273.5521676 28.0543278 23.9453268 443.4580780 26.9993527 728.9399978 18.7638775 301.0758970 25.3739780 25.0186104 332.9325937 3.5626571 67.8801268 2.0528618 303.3227320 6.9689095 3.5529388 54.2456364 4.1215053 78.6792379 5.150162 70.7753041 4.4225772 4.2561245 446.1703598 20.8869507 
write.table(de3,file="de3b.tsv",quote=FALSE,sep="\t")

ss4 %>% kbl() %>% kable_paper("hover", full_width = F)
UDI fraction Group RatID label nreads mtfrac trt
1 1 wholetissue T 1 wholetissue T 22353545 0.2350623 0
10 10 wholetissue T 10 wholetissue T 32594141 0.3531781 0
13 13 wholetissue T 13 wholetissue T 29989340 0.2972082 0
14 14 wholetissue T 14 wholetissue T 27775661 0.3282823 0
17 17 wholetissue T 17 wholetissue T 30401193 0.3348082 0
18 18 wholetissue T 18 wholetissue T 33702596 0.3547981 0
2 2 wholetissue T 2 wholetissue T 19477976 0.3181207 0
21 21 wholetissue T 21 wholetissue T 39401323 0.3400537 0
22 22 wholetissue T 22 wholetissue T 27353327 0.3221241 0
25 25 wholetissue T 25 wholetissue T 34339726 0.2938285 0
26 26 wholetissue T 26 wholetissue T 28163757 0.2587561 0
27 27 wholetissue T 27 wholetissue T 30543678 0.3380834 0
28 28 wholetissue T 28 wholetissue T 30906624 0.2776614 0
5 5 wholetissue T 5 wholetissue T 30755504 0.3586038 0
6 6 wholetissue T 6 wholetissue T 26148398 0.3245408 0
9 9 wholetissue T 9 wholetissue T 24737507 0.3535222 0
33 33 mito T 1 mito T 5474018 0.8850587 1
34 34 mito T 2 mito T 9988027 0.9542655 1
37 37 mito T 5 mito T 1353729 0.4835828 1
46 46 mito T 14 mito T 1086442 0.4345745 1
53 53 mito T 21 mito T 1628001 0.5591495 1
57 57 mito T 25 mito T 1841975 0.5667390 1
58 58 mito T 26 mito T 3802593 0.7906639 1
59 59 mito T 27 mito T 7532868 0.9057456 1
60 60 mito T 28 mito T 6001091 0.8726930 1 
de4 <- run_de(ss4,x4)
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 420 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
##  chr [1:661] "ENSRNOG00000006049 Rfx1 protein_coding" ...
##  chr [1:77] "ENSRNOG00000001052 Slc25a30 protein_coding" ...

as.data.frame(de4[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)
baseMean log2FoldChange lfcSE stat pvalue padj 1 10 13 14 17 18 2 21 22 25 26 27 28 33 34 37 46 5 53 57 58 59 6 60 9
ENSRNOG00000000001 Arsj protein_coding 35.629246 1.3761157 0.5839760 2.3564594 0.0184501 0.1006368 0.8052414 1.1362120 8.1523318 0.5340835 36.1963140 0.8947127 0.7811458 4.2076551 0.5340835 14.7740057 0.5368276 0.9586789 4.2076551 0.4153983 8.8644034 1.2525977 0.4970928 4.9965905 0.9198105 15.5127060 1.1631265 0.4970928 4.7336120 0.7121113 9.6031037
ENSRNOG00000000007 Gad1 protein_coding 57.291868 1.9204448 0.8808192 2.1802941 0.0292357 0.1181549 0.2454429 0.1637000 1.5930188 0.0513400 8.2839213 0.2454429 0.0982200 2.7877828 0.2567002 11.0452284 0.1534018 0.3601400 1.3275156 1.7455612 19.3291496 1.1351733 1.2441200 1.0620125 1.6942212 53.3852705 1.1658537 1.1131600 2.5222797 2.5156618 5.5226142
ENSRNOG00000000008 Alx4 protein_coding 30.766047 0.7797296 0.4527809 1.7220903 0.0850532 0.1942579 0.2000711 1.4883541 0.9560815 1.1167138 2.3735589 1.3338073 0.6794660 2.4475687 0.2791784 1.1380077 1.4671880 2.0060425 1.1090546 0.3553180 0.3251451 0.5668681 0.2911997 0.9560815 0.3299382 0.6502901 0.4001422 0.3882663 1.0708113 0.5583569 0.8128626
ENSRNOG00000000009 Tmco5b protein_coding 24.107869 1.7255317 1.7108656 1.0085723 0.3131798 0.4214014 0.0360027 0.1826812 0.0000000 0.0365586 0.0000000 0.0000000 0.0000000 0.3332727 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.6580552 4.2997517 0.8280631 2.3748552 0.0000000 0.8042897 17.1990066 0.5400412 2.7402175 0.0000000 0.3655862 0.0000000
ENSRNOG00000000010 Cbln1 protein_coding 34.054233 1.7828863 1.5589297 1.1436605 0.2527645 0.3745620 0.0000000 0.0000000 0.0000000 0.0291208 0.5428955 0.0000000 0.4004795 0.0404244 0.1747247 0.0000000 0.0000000 0.0000000 0.0808489 0.6697782 5.9718509 0.9868034 2.5029968 0.0000000 0.6988990 9.2292241 0.8881230 2.7032366 0.0000000 0.7280198 0.0000000
ENSRNOG00000000012 Tcf15 protein_coding 180.539920 -0.7432341 0.4705600 -1.5794672 0.1142289 0.2309516 8.2486227 500.8387942 31.4044148 15.8359554 96.7760683 14.3015689 237.8614922 40.5304841 21.2329626 162.7836584 19.2270055 640.4531483 33.9096103 0.6746259 2.8927629 0.1483565 7.3870029 37.2200472 0.3550663 3.1557414 0.2373704 3.6935014 28.8992193 0.3550663 115.1845596
ENSRNOG00000000017 Steap1 protein_coding 12.487775 0.9320222 1.7524382 0.5318431 0.5948346 0.6349797 0.0000000 0.0000000 0.0306804 0.0000000 0.3982547 0.0000000 0.9204357 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.3601400 0.5310063 1.0268007 5.5226142 0.0000000 0.2291800 1.4602672 0.1540201 5.5226142 0.0000000 0.2946600 0.0000000
ENSRNOG00000000021 AABR07061902.1 protein_coding 11.214462 1.8176431 1.9513598 0.9314751 0.3516079 0.4482125 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0253799 0.0000000 0.0000000 0.0970666 0.1529730 0.3806979 0.2601160 0.0000000 0.0970666 0.4971624 0.3806979 0.2601160 0.0000000 0.0647110 0.0000000
ENSRNOG00000000024 Hebp1 protein_coding 30.006590 0.0991046 0.1734657 0.5713212 0.5677819 0.6120443 1.4623450 56.5110218 3.0602332 9.4994207 12.6643639 2.2300761 22.1130085 2.7362085 13.7010876 9.4982729 2.3763106 47.9115185 2.6282003 0.9134058 0.9998182 0.2193517 5.5282521 1.9801509 1.0960870 1.8330000 0.4021449 6.1425024 2.4481866 2.0094928 6.4988183
ENSRNOG00000000033 Tmcc2 protein_coding 350.908906 -0.8779014 0.2868408 -3.0605879 0.0022090 0.0555430 30.8389182 546.6958021 37.5972088 99.1186748 42.7690632 34.7411042 346.3673503 40.3273648 108.6300628 54.8155479 35.5273656 726.9371191 43.3864553 3.1037161 0.9297622 1.3977980 21.1729258 30.9198392 3.9046751 1.7382512 1.2230732 16.8297615 30.8211589 2.2026372 25.8716450
ENSRNOG00000000034 Nuak2 protein_coding 68.826907 1.1665124 0.4264896 2.7351486 0.0062352 0.0729163 2.1659042 21.8272111 2.5220609 39.1511152 2.7736093 1.3492518 10.2561594 1.7802783 54.6638212 3.8920001 2.4854639 29.9795429 1.2165235 25.1158097 1.0736552 0.9586789 8.6782887 2.7297600 18.4675072 1.7894253 1.6688114 10.5191379 1.6022505 28.0706109 2.4157242
ENSRNOG00000000035 Tmed11 protein_coding 6.673576 1.3150727 2.0391199 0.6449217 0.5189779 0.5708222 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 4.6021785 0.1840822 0.0982200 0.3982547 0.0000000 3.6817428 0.2147625 0.0654800 0.5310063 0.0000000 4.6021785 0.0000000
ENSRNOG00000000036 Klhdc8a protein_coding 44.759178 2.1376054 1.1836045 1.8060132 0.0709163 0.1762416 0.0323555 0.2677028 0.0253799 0.1950870 0.2667615 0.0000000 0.0000000 0.0507597 0.0975435 0.1667259 0.0647110 0.0764865 0.0253799 1.0079497 0.4668325 0.7441770 0.9943248 0.0000000 1.4306382 1.3671525 1.1971544 1.2237843 0.0000000 0.9104062 0.0666904
ENSRNOG00000000040 RGD1304622 protein_coding 106.644402 0.8961277 0.3311529 2.7060842 0.0068082 0.0752479 21.0083343 28.1615460 5.7031454 74.9385289 4.5363457 27.2194940 13.8308184 5.4472350 70.6387772 4.5003429 24.6619576 24.8288186 4.9354143 36.8550142 0.7920604 8.2206525 6.4988183 6.0321730 15.9705062 2.6642030 9.8647831 10.3314547 5.2644419 30.7125118 3.6362771
ENSRNOG00000000041 Slc26a1 protein_coding 42.451276 2.0868257 1.4331467 1.4561144 0.1453610 0.2671025 1.1013186 0.0404244 0.0000000 0.5428955 0.0000000 0.1001199 0.3233956 0.0000000 1.6286866 0.0000000 0.0000000 0.0404244 0.0000000 17.3726571 0.2302541 3.1037161 1.0914600 0.0000000 19.0013437 1.2170575 4.3051546 1.0106111 0.0291208 16.2868660 0.0000000
ENSRNOG00000000042 Xpr1 protein_coding 194.044688 -0.2099177 0.1078470 -1.9464406 0.0516018 0.1502819 342.0182326 22.9493801 17.8953397 118.6032794 14.9543376 386.3402498 14.6285522 23.5053867 132.8041155 16.6752733 378.2145466 22.7704375 19.4931379 6.0485043 0.7417826 31.7641123 1.4762759 14.5577169 11.5710516 1.2461948 30.2867117 1.5210115 14.9482897 19.9863619 10.4146280
ENSRNOG00000000043 Idua protein_coding 150.095911 -0.2251701 0.1069270 -2.1058309 0.0352190 0.1273399 274.2898378 11.5664960 14.1764197 44.6045251 17.7123126 370.9355861 8.7439028 13.8817598 50.7110970 23.0003364 336.8794653 14.2970480 14.6347797 3.9825469 1.2321609 31.2948137 1.1351733 13.0632598 3.4515406 1.6428812 28.5335066 0.7670090 11.3935198 3.7170438 14.1698501
ENSRNOG00000000044 Tmem175 protein_coding 117.276824 -0.1311654 0.1088068 -1.2054886 0.2280147 0.3524287 7.3807927 10.7038034 9.2213242 10.2491939 6.8271819 9.1040615 5.7687165 13.3304757 10.2109506 9.7458656 8.1936554 11.7375041 10.9361670 0.7648652 0.6344965 0.7478336 0.7002488 8.6712803 0.9178383 0.9898145 0.6828046 1.0337006 8.4124361 0.7266220 6.2180653
ENSRNOG00000000047 Cd82 protein_coding 108.242117 0.8122696 0.2812947 2.8876105 0.0038818 0.0640084 76.7812796 5.8684472 27.2194940 20.6629095 7.3848421 92.7517858 3.8522936 41.8339874 24.6621823 5.4472350 84.1522825 6.5164966 29.9597115 7.8319092 0.5118207 33.1695128 1.5121152 34.1613784 11.1646366 2.3763106 29.4840114 1.5481180 23.9312330 5.3323637 4.2408004
ENSRNOG00000000048 Gak protein_coding 316.310469 0.1641226 0.0831985 1.9726633 0.0485339 0.1458256 276.3338265 23.3543467 76.1912238 22.0313227 18.2878570 448.4317105 12.3679357 82.3986559 30.1162118 23.8499282 386.5416197 24.6700845 78.3938610 2.7084378 1.6598851 47.7748069 3.0590905 66.7799556 3.3956534 3.5527366 66.7761506 2.3354347 61.6738421 3.8807468 14.5895165 
write.table(de4,file="de4b.tsv",quote=FALSE,sep="\t")

ss5 %>% kbl() %>% kable_paper("hover", full_width = F)
UDI fraction Group RatID label nreads mtfrac trt
34 34 mito T 2 mito T 9988027 0.9542655 1
47 47 mito S 15 mito S 8344505 0.9222841 0
59 59 mito T 27 mito T 7532868 0.9057456 1
60 60 mito T 28 mito T 6001091 0.8726930 1 
de5 <- run_de_cov(ss5,x5)
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.
## final dispersion estimates
## fitting model and testing
##  chr [1:43] "ENSRNOG00000058403 AABR07060960.1 snoRNA" ...
##  chr(0)

as.data.frame(de5[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)
baseMean log2FoldChange lfcSE stat pvalue padj 34 47 59 60
ENSRNOG00000000001 Arsj protein_coding 13.879255 0.9925647 1.0906845 0.9100383 0.3628023 0.9999545 1.2014385 0.8009590 1.4016782 2.4028770
ENSRNOG00000000007 Gad1 protein_coding 34.371879 -0.3771957 0.6771853 -0.5570051 0.5775240 0.9999545 2.5166262 4.9134131 4.0745377 5.8721278
ENSRNOG00000000008 Alx4 protein_coding 17.989689 -1.1840777 0.9112336 -1.2994228 0.1937989 0.9999545 1.3275156 4.1152984 1.5930188 2.9205344
ENSRNOG00000000009 Tmco5b protein_coding 12.807227 -0.9773889 1.0669809 -0.9160322 0.3596500 0.9999545 1.1664546 3.4993637 2.4995455 1.6663637
ENSRNOG00000000010 Cbln1 protein_coding 19.779707 -0.1423100 0.8830251 -0.1611619 0.8719659 0.9999545 1.1013186 2.1025173 2.7032366 2.5029968
ENSRNOG00000000012 Tcf15 protein_coding 8.299133 0.6973686 1.4115390 0.4940484 0.6212720 0.9999545 1.3182328 0.7190361 0.5991967 1.1983934
ENSRNOG00000000017 Steap1 protein_coding 6.422884 -0.4558564 1.5567867 -0.2928188 0.7696606 0.9999545 0.5310063 1.0620125 0.7965094 1.1947641
ENSRNOG00000000021 AABR07061902.1 protein_coding 4.255812 0.7945281 2.1190525 0.3749450 0.7077014 0.9999545 0.6665455 0.4999091 1.3330909 0.3332727
ENSRNOG00000000024 Hebp1 protein_coding 8.136685 0.3088650 1.3782218 0.2241040 0.8226764 0.9999545 0.6007192 0.7008391 1.0011987 1.1013186
ENSRNOG00000000033 Tmcc2 protein_coding 25.874495 -0.0337981 0.7719970 -0.0437800 0.9650797 0.9999545 2.7563049 3.3555016 3.7150196 2.6364656
ENSRNOG00000000034 Nuak2 protein_coding 31.605287 0.2181646 0.7102859 0.3071504 0.7587289 0.9999545 3.1860375 3.8497953 5.3100625 5.0445594
ENSRNOG00000000035 Tmed11 protein_coding 4.650448 0.9026108 1.9776182 0.4564131 0.6480929 0.9999545 0.9998182 0.4999091 0.6665455 0.8331818
ENSRNOG00000000036 Klhdc8a protein_coding 24.455051 -0.3306512 0.7980813 -0.4143077 0.6786488 0.9999545 1.4016782 2.9034763 3.2038360 2.8033565
ENSRNOG00000000040 RGD1304622 protein_coding 45.725955 -0.5753018 0.6168786 -0.9326013 0.3510259 0.9999545 2.6364656 7.1903606 7.4300393 5.9919672
ENSRNOG00000000041 Slc26a1 protein_coding 23.123348 -1.2829039 0.8238276 -1.5572480 0.1194116 0.9999545 0.9292609 5.0445594 3.3187891 3.9825469
ENSRNOG00000000042 Xpr1 protein_coding 37.892972 0.4437841 0.6900713 0.6430989 0.5201599 0.9999545 4.1659092 4.6658183 5.6656365 12.6643639
ENSRNOG00000000043 Idua protein_coding 26.080652 -0.0367088 0.7679198 -0.0478030 0.9618733 0.9999545 2.4028770 2.8033565 2.5029968 2.8033565
ENSRNOG00000000044 Tmem175 protein_coding 26.208470 -0.0663783 0.7689298 -0.0863255 0.9312076 0.9999545 2.9959836 3.4753410 3.7150196 2.2769475
ENSRNOG00000000047 Cd82 protein_coding 31.338805 -0.7236387 0.7314654 -0.9893001 0.3225164 0.9999545 1.8585219 5.8410688 5.7083172 4.2480500
ENSRNOG00000000048 Gak protein_coding 77.552594 -0.3757874 0.4839096 -0.7765653 0.4374153 0.9999545 9.4982729 15.9970912 11.8311820 15.9970912 
write.table(de5,file="de5b.tsv",quote=FALSE,sep="\t")

Top expressed genes in low medium and high mtfrac samples

# high
high <- x[which(colnames(x) %in% subset(mito,mtfrac>0.6)$RatID )] 
high <- high/colSums(high)*1000000 
high$mean <- rowMeans(high)
head(high[order(-high$mean),],20) %>% kbl() %>% kable_paper("hover", full_width = F)
1 15 2 20 24 26 27 28 32 7 mean
ENSRNOG00000043866 AY172581.24 Mt_rRNA 94682.78429 150009.221 104054.600 215428.646 140492.986 146753.133 212385.671 170947.895 135245.777 139515.369 150951.608
ENSRNOG00000030478 AY172581.9 Mt_rRNA 65586.97579 90635.554 63821.561 125909.740 88687.784 92318.268 122501.280 104777.474 85580.230 89331.918 92915.079
ENSRNOG00000069271 Ttn protein_coding 59601.46129 68025.573 44342.954 80628.887 86504.771 80265.162 55496.241 68918.833 82118.699 57774.340 68367.692
ENSRNOG00000049695 Myh4 protein_coding 75056.90932 50916.040 29081.070 71544.849 75531.599 62216.396 28406.546 70644.768 30227.510 60595.758 55422.144
ENSRNOG00000017786 Acta1 protein_coding 28957.84068 31772.665 24726.001 29042.431 28878.279 34457.537 34758.970 41810.807 18510.875 29637.872 30255.328
ENSRNOG00000034234 Mt-co1 protein_coding 18182.96382 26811.615 30763.300 24765.916 20807.909 22750.709 37794.658 28969.747 18940.098 29843.296 25963.021
ENSRNOG00000006783 Neb protein_coding 15534.85751 15830.975 13376.226 16532.799 16618.670 18201.406 15552.847 17875.586 15404.206 13666.747 15859.432
ENSRNOG00000047124 AABR07005775.1 protein_coding 15960.38518 13545.029 7877.141 18413.835 20063.934 18357.034 14099.716 20942.992 7469.084 14559.349 15128.850
ENSRNOG00000065740 Myh2 protein_coding 1865.85095 9661.681 9632.572 13552.589 9840.961 16110.093 27439.915 15396.210 8463.662 9408.113 12137.165
ENSRNOG00000016837 Ckm protein_coding 10450.89374 10599.797 7271.796 10168.365 11130.590 11733.622 10861.783 15001.151 6294.294 11137.252 10464.954
ENSRNOG00000018184 Tpm1 protein_coding 14317.86342 9205.217 6630.771 12278.945 12191.873 10924.595 5136.314 12453.399 5750.341 12296.195 10118.551
ENSRNOG00000020332 Tnnt3 protein_coding 10081.49123 8755.197 6353.795 9793.710 9970.412 9351.069 7095.138 10774.264 5529.059 10335.599 8803.974
ENSRNOG00000016983 Myh7 protein_coding 34.96676 8904.357 14843.149 5629.169 4111.660 9149.747 20893.548 8260.921 11409.318 4151.608 8738.844
ENSRNOG00000029707 Mt-nd4 protein_coding 7365.53210 9067.599 10770.911 7901.960 6369.842 7000.345 11290.718 9245.498 6001.541 10110.423 8512.437
ENSRNOG00000017645 Mylpf protein_coding 10032.76726 8433.063 7113.367 9993.543 8025.269 6852.201 4942.351 9460.634 5316.877 11071.684 8124.176
ENSRNOG00000031766 Mt-cyb protein_coding 6600.78850 8192.381 10356.158 7433.524 6370.426 7227.548 10273.449 8982.966 5473.053 8597.043 7950.734
ENSRNOG00000020557 Ryr1 protein_coding 6877.07475 6956.928 3852.678 6634.966 8412.491 8784.065 6074.160 9280.184 4405.520 6374.873 6765.294
ENSRNOG00000023803 Cmya5 protein_coding 5672.94939 5545.919 3160.947 6048.339 6897.936 6174.794 4906.315 6867.519 3937.832 5349.196 5456.175
ENSRNOG00000031979 Mt-atp6 protein_coding 4067.10247 6027.343 5879.513 5379.106 4245.074 4709.750 6912.855 6037.608 4124.672 6580.841 5396.386
ENSRNOG00000058068 Obscn protein_coding 4199.20835 4938.982 2724.337 4224.806 5894.355 6190.996 5059.351 6229.376 3621.846 4577.284 4766.054 
# med
med <- x[which(colnames(x) %in% subset(mito,mtfrac<0.6 & mtfrac > 0.3)$RatID )]  
med <- med/colSums(med)*1000000
med$mean <- rowMeans(med)
head(med[order(-med$mean),],20) %>% kbl() %>% kable_paper("hover", full_width = F)
11 14 16 18 21 22 25 29 3 30 31 4 5 mean
ENSRNOG00000030478 AY172581.9 Mt_rRNA 145096.173 121452.685 111027.988 154442.242 210917.201 111294.598 124266.392 116839.486 105031.494 119309.763 99219.105 88165.787 155818.030 127913.919
ENSRNOG00000043866 AY172581.24 Mt_rRNA 137369.699 83494.709 108521.774 193082.418 217162.951 121757.547 128802.422 106172.671 108134.183 75711.961 105691.354 87757.079 136185.734 123834.192
ENSRNOG00000069271 Ttn protein_coding 62241.797 74560.825 97344.143 76348.934 111391.264 83840.653 118473.654 50725.603 48863.658 87257.552 68894.933 75841.316 45199.100 76998.726
ENSRNOG00000049695 Myh4 protein_coding 32275.784 34743.415 44333.784 31826.824 56469.184 35765.457 48358.589 23968.058 40562.014 51707.271 20354.266 62303.468 22871.902 38887.694
ENSRNOG00000034234 Mt-co1 protein_coding 43692.196 25548.025 33131.633 35920.766 46831.735 25855.465 35938.552 28001.687 24628.768 21060.797 25945.093 23033.283 52121.574 32439.198
ENSRNOG00000017786 Acta1 protein_coding 28723.981 18198.348 24523.413 25639.578 32919.327 18664.128 27014.942 19980.229 24454.731 21644.116 16921.219 24872.650 25989.488 23811.242
ENSRNOG00000006783 Neb protein_coding 17802.368 18803.206 17518.458 26266.144 19286.438 13079.319 24375.362 14094.360 15811.509 13696.215 14509.970 21647.009 18834.666 18132.694
ENSRNOG00000047124 AABR07005775.1 protein_coding 9418.063 14146.395 9002.388 13332.248 20909.696 18060.699 20815.749 9165.848 10825.932 15429.228 8380.789 11951.097 10291.338 13209.959
ENSRNOG00000065740 Myh2 protein_coding 11523.242 12056.458 7986.570 14721.398 11664.391 10004.749 12208.757 10018.991 5660.310 6798.883 11087.765 3941.840 21450.500 10701.835
ENSRNOG00000016837 Ckm protein_coding 8892.157 10617.143 11931.659 6648.355 11105.730 8048.955 12757.294 8396.200 10039.359 9922.126 8326.240 8720.337 7905.593 9485.473
ENSRNOG00000016983 Myh7 protein_coding 16133.752 9364.955 9398.950 11555.330 9617.398 3055.479 9370.523 10019.780 5228.963 9284.797 7479.456 3045.475 11546.356 8853.939
ENSRNOG00000029707 Mt-nd4 protein_coding 7812.854 8518.652 9262.162 5074.144 8492.759 5536.513 8397.306 8905.864 5577.157 5504.283 7677.849 4618.528 8675.682 7234.904
ENSRNOG00000031979 Mt-atp6 protein_coding 7932.953 8698.326 8622.431 7462.200 8533.314 4825.153 8278.297 5820.518 4904.993 5727.293 7352.659 6136.751 7618.196 7070.237
ENSRNOG00000020332 Tnnt3 protein_coding 5780.097 6653.988 8099.445 4650.377 8504.424 7261.610 8837.684 5214.101 6529.541 6099.396 5588.925 7070.120 5083.436 6567.165
ENSRNOG00000031766 Mt-cyb protein_coding 6906.223 6746.857 7068.046 5724.915 6474.664 5385.506 7642.408 6726.645 3338.670 5271.306 6280.022 5463.026 10395.043 6417.179
ENSRNOG00000018184 Tpm1 protein_coding 5637.996 5927.120 4035.550 6442.259 9457.635 8373.212 7734.854 4339.390 7121.555 6489.390 3846.682 9099.744 4145.448 6357.756
ENSRNOG00000017645 Mylpf protein_coding 5376.001 3801.661 3891.718 4756.865 7401.586 5135.412 3662.440 5090.464 8296.443 5904.575 5832.088 9177.553 4900.879 5632.899
ENSRNOG00000015155 Tnnc2 protein_coding 5849.130 5584.517 6060.320 3263.799 6457.524 4545.012 5400.491 4750.574 6748.189 4475.004 5152.782 6297.400 5379.675 5381.878
ENSRNOG00000058068 Obscn protein_coding 3636.692 4813.843 6019.366 6218.383 8180.442 3438.463 7367.493 4489.791 4084.208 6387.267 3902.310 4202.402 4465.391 5169.696
ENSRNOG00000020557 Ryr1 protein_coding 4119.409 4815.463 6062.903 4455.103 5844.179 5293.179 8230.124 4460.859 2878.837 6427.246 4159.975 6328.305 3863.552 5149.164 
# low
low <- x[which(colnames(x) %in% subset(mito,mtfrac<0.3)$RatID )]
low <- low/colSums(low)*1000000
low$mean <- rowMeans(low)
head(low[order(-low$mean),],20) %>% kbl() %>% kable_paper("hover", full_width = F)
10 12 13 17 19 23 6 8 9 mean
ENSRNOG00000030478 AY172581.9 Mt_rRNA 168395.425 93179.371 97805.063 135794.378 140435.707 183670.564 108069.801 85452.902 100148.372 123661.287
ENSRNOG00000043866 AY172581.24 Mt_rRNA 138016.929 112963.823 109505.278 179749.662 128139.281 173566.884 102244.649 59991.917 99581.427 122639.983
ENSRNOG00000069271 Ttn protein_coding 58854.313 98810.827 68381.066 91817.377 85235.585 69612.296 57264.535 98383.944 52659.058 75668.778
ENSRNOG00000049695 Myh4 protein_coding 23359.215 66938.152 36668.016 28555.439 37161.248 23037.404 31270.137 53197.836 37690.364 37541.979
ENSRNOG00000034234 Mt-co1 protein_coding 42774.320 27999.239 38469.661 26058.793 35413.720 33640.292 36566.296 13226.186 31541.344 31743.317
ENSRNOG00000017786 Acta1 protein_coding 25814.969 27509.480 30732.246 20018.413 25884.853 23444.508 24997.507 18817.306 18810.129 24003.268
ENSRNOG00000006783 Neb protein_coding 18021.521 14724.425 16268.982 20283.724 22314.598 18791.010 20514.415 16972.623 16701.764 18288.118
ENSRNOG00000047124 AABR07005775.1 protein_coding 11201.986 15961.438 17469.895 11362.449 17040.763 9746.084 9360.639 11621.366 8954.201 12524.314
ENSRNOG00000065740 Myh2 protein_coding 17891.981 5017.008 11052.845 13087.642 12983.468 18003.206 12255.988 2936.636 9359.731 11398.723
ENSRNOG00000016983 Myh7 protein_coding 20103.144 5635.703 9495.943 12735.553 12169.294 15672.297 8768.232 2090.322 8426.910 10566.377
ENSRNOG00000016837 Ckm protein_coding 9073.717 11993.347 10920.554 8321.963 7644.216 10535.026 8006.034 11584.595 8054.574 9570.447
ENSRNOG00000029707 Mt-nd4 protein_coding 9423.529 6359.732 6864.719 5094.863 5997.423 7026.463 8753.084 3412.794 10373.438 7034.005
ENSRNOG00000031979 Mt-atp6 protein_coding 11286.465 6266.632 6736.333 4950.187 5523.764 7549.296 6562.665 4385.767 7356.257 6735.263
ENSRNOG00000020332 Tnnt3 protein_coding 5970.155 10596.406 7103.965 6585.344 5921.730 4776.503 5133.058 6762.151 6591.723 6604.559
ENSRNOG00000031766 Mt-cyb protein_coding 8744.574 7827.157 6741.676 6704.192 5352.134 5742.873 5588.333 3237.582 7404.278 6371.422
ENSRNOG00000018184 Tpm1 protein_coding 4551.704 7557.616 6216.840 6004.738 6549.235 4314.166 6638.646 6883.283 5070.842 5976.341
ENSRNOG00000015155 Tnnc2 protein_coding 7425.328 8285.004 7999.755 3587.545 3605.801 3738.595 4173.510 5566.039 5125.421 5500.778
ENSRNOG00000058068 Obscn protein_coding 4969.155 7213.750 6094.532 6803.979 5212.493 6293.844 4243.764 4779.110 3443.980 5450.512
ENSRNOG00000020557 Ryr1 protein_coding 5336.633 6114.790 5124.541 4207.971 5089.772 5326.215 4083.751 8176.809 3519.566 5220.005
ENSRNOG00000017645 Mylpf protein_coding 4537.027 9098.033 6005.667 4879.309 3714.394 4580.656 4583.242 5005.317 4277.382 5186.781

Pathway analysis with mitch

Here we are doing a gene set analysis with my R package called mitch. I’m using gene sets downloaded from Reactome 5th Dec 2020.

We lost 46% of genes after converting from rat to human.

There were 259 differentially regulated gene sets (FDR<0.05). Of these 81 were downregulated and 178 were upregulated.

#download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip",
#    destfile="ReactomePathways.gmt.zip")
#unzip("ReactomePathways.gmt.zip")
genesets<-gmt_import("ReactomePathways.gmt")

# i need to get some data from biomart to link rat and human gene IDs

mart <- read.table("mart_export.txt",header=TRUE,sep="\t")
gt <- mart[,c("Gene.stable.ID","Human.gene.name")]

rownames(de1) <- sapply( strsplit(rownames(de1)," ") , "[[", 1)

m <- mitch_import(as.data.frame(de1),"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 = 30560
## Note: no. genes in output = 16477
## Note: estimated proportion of input genes in output = 0.539
res<-mitch_calc(m,genesets,priority="significance")
## Note: When prioritising by significance (ie: small
##             p-values), large effect sizes might be missed.
nrow(subset(res$enrichment_result,p.adjustANOVA<0.05))
## [1] 79
head(res$enrichment_result,20) %>% kbl() %>% kable_paper("hover", full_width = F)
set setSize pANOVA s.dist p.adjustANOVA
648 Major pathway of rRNA processing in the nucleolus and cytosol 148 0 -0.3703409 0e+00
1073 SRP-dependent cotranslational protein targeting to membrane 83 0 -0.4837848 0e+00
1396 rRNA processing in the nucleus and cytosol 155 0 -0.3544727 0e+00
1040 Response of EIF2AK4 (GCN2) to amino acid deficiency 72 0 -0.4845863 0e+00
833 Peptide chain elongation 60 0 -0.5268776 0e+00
365 Eukaryotic Translation Termination 63 0 -0.5137803 0e+00
1353 Viral mRNA Translation 60 0 -0.5255518 0e+00
1093 Selenocysteine synthesis 63 0 -0.5121886 0e+00
655 Metabolism 1845 0 0.0996703 0e+00
363 Eukaryotic Translation Elongation 63 0 -0.5089374 0e+00
149 Cap-dependent Translation Initiation 87 0 -0.4225923 0e+00
364 Eukaryotic Translation Initiation 87 0 -0.4225923 0e+00
765 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) 65 0 -0.4809773 0e+00
456 GTP hydrolysis and joining of the 60S ribosomal subunit 81 0 -0.4309866 0e+00
415 Formation of a pool of free 40S subunits 71 0 -0.4562587 0e+00
610 L13a-mediated translational silencing of Ceruloplasmin expression 79 0 -0.4307186 0e+00
1265 The citric acid (TCA) cycle and respiratory electron transport 150 0 0.3094620 0e+00
764 Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC) 83 0 -0.4089984 0e+00
766 Nonsense-Mediated Decay (NMD) 83 0 -0.4089984 0e+00
1394 rRNA processing 176 0 -0.2645996 1e-07 
mitch_barplot <- function(res){
  sig <- head(subset(res$enrichment_result,p.adjustANOVA<0.05),30)
  sig <- sig[order(sig$s.dist),]
  par(mar=c(3,25,1,1)); barplot(sig$s.dist,horiz=TRUE,las=2,cex.names = 0.6,cex.axis = 0.6,
    names.arg=sig$set,main="Enrichment score") ;grid()
}

mitch_barplot(res)

nrow(subset(res$enrichment_result,p.adjustANOVA<0.05&s.dist<0))
## [1] 51
nrow(subset(res$enrichment_result,p.adjustANOVA<0.05&s.dist>0))
## [1] 28
unlink("skeletal_mitch1.html")
mitch_report(res,outfile="skeletal_mitch1.html")
## Dataset saved as " /tmp/RtmpwUhoHt/skeletal_mitch1.rds ".
## 
## 
## processing file: mitch.Rmd
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |..                                                                    |   3%
##    inline R code fragments
## 
## 
  |                                                                            
  |....                                                                  |   6%
## label: checklibraries (with options) 
## List of 1
##  $ echo: logi FALSE
## 
## 
  |                                                                            
  |......                                                                |   9%
##   ordinary text without R code
## 
## 
  |                                                                            
  |........                                                              |  12%
## label: peek (with options) 
## List of 1
##  $ echo: logi FALSE
## 
## 
  |                                                                            
  |..........                                                            |  15%
##   ordinary text without R code
## 
## 
  |                                                                            
  |............                                                          |  18%
## label: metrics (with options) 
## List of 1
##  $ echo: logi FALSE
## 
## 
  |                                                                            
  |..............                                                        |  21%
##   ordinary text without R code
## 
## 
  |                                                                            
  |................                                                      |  24%
## label: scatterplot (with options) 
## List of 5
##  $ echo      : logi FALSE
##  $ fig.height: num 6
##  $ fig.width : num 6.5
##  $ message   : logi FALSE
##  $ warning   : logi FALSE
## 
  |                                                                            
  |...................                                                   |  26%
##   ordinary text without R code
## 
## 
  |                                                                            
  |.....................                                                 |  29%
## label: contourplot (with options) 
## List of 5
##  $ echo      : logi FALSE
##  $ fig.height: num 6
##  $ fig.width : num 6.5
##  $ warning   : logi FALSE
##  $ message   : logi FALSE
## Contour plot does not apply to unidimensional analysis.
## 
  |                                                                            
  |.......................                                               |  32%
##   ordinary text without R code
## 
## 
  |                                                                            
  |.........................                                             |  35%
## label: input_geneset_metrics1 (with options) 
## List of 2
##  $ results: chr "asis"
##  $ echo   : logi FALSE
## 
## 
  |                                                                            
  |...........................                                           |  38%
##   ordinary text without R code
## 
## 
  |                                                                            
  |.............................                                         |  41%
## label: input_geneset_metrics2 (with options) 
## List of 5
##  $ results   : chr "asis"
##  $ echo      : logi FALSE
##  $ fig.height: num 7
##  $ fig.width : num 7
##  $ fig.show  : chr "all"
## 
  |                                                                            
  |...............................                                       |  44%
##   ordinary text without R code
## 
## 
  |                                                                            
  |.................................                                     |  47%
## label: input_geneset_metrics3 (with options) 
## List of 5
##  $ results   : chr "asis"
##  $ echo      : logi FALSE
##  $ message   : logi FALSE
##  $ fig.height: num 7
##  $ fig.width : num 7
## 
  |                                                                            
  |...................................                                   |  50%
##   ordinary text without R code
## 
## 
  |                                                                            
  |.....................................                                 |  53%
## label: echart1d (with options) 
## List of 6
##  $ results   : chr "asis"
##  $ echo      : logi FALSE
##  $ fig.height: num 7
##  $ fig.width : num 7
##  $ fig.show  : chr "all"
##  $ message   : logi FALSE
## 
## 
  |                                                                            
  |.......................................                               |  56%
## label: echart2d (with options) 
## List of 6
##  $ results   : chr "asis"
##  $ echo      : logi FALSE
##  $ fig.height: num 7
##  $ fig.width : num 7
##  $ fig.show  : chr "all"
##  $ message   : logi FALSE
## 
## 
  |                                                                            
  |.........................................                             |  59%
##   ordinary text without R code
## 
## 
  |                                                                            
  |...........................................                           |  62%
## label: heatmap (with options) 
## List of 6
##  $ results   : chr "asis"
##  $ echo      : logi FALSE
##  $ fig.height: num 10
##  $ fig.width : num 7
##  $ fig.show  : chr "all"
##  $ message   : logi FALSE
## 
## 
  |                                                                            
  |.............................................                         |  65%
##   ordinary text without R code
## 
## 
  |                                                                            
  |...............................................                       |  68%
## label: effectsize (with options) 
## List of 6
##  $ results   : chr "asis"
##  $ echo      : logi FALSE
##  $ fig.height: num 7
##  $ fig.width : num 7
##  $ fig.show  : chr "all"
##  $ message   : logi FALSE
## 
## 
  |                                                                            
  |.................................................                     |  71%
##   ordinary text without R code
## 
## 
  |                                                                            
  |...................................................                   |  74%
## label: results_table (with options) 
## List of 2
##  $ results: chr "asis"
##  $ echo   : logi FALSE
## 
## 
  |                                                                            
  |......................................................                |  76%
##   ordinary text without R code
## 
## 
  |                                                                            
  |........................................................              |  79%
## label: results_table_complete (with options) 
## List of 2
##  $ results: chr "asis"
##  $ echo   : logi FALSE
## 
## 
  |                                                                            
  |..........................................................            |  82%
##   ordinary text without R code
## 
## 
  |                                                                            
  |............................................................          |  85%
## label: detailed_geneset_reports1d (with options) 
## List of 7
##  $ results   : chr "asis"
##  $ echo      : logi FALSE
##  $ fig.height: num 6
##  $ fig.width : num 6
##  $ out.width : chr "80%"
##  $ comment   : logi NA
##  $ message   : logi FALSE
## 
  |                                                                            
  |..............................................................        |  88%
##   ordinary text without R code
## 
## 
  |                                                                            
  |................................................................      |  91%
## label: detailed_geneset_reports2d (with options) 
## List of 7
##  $ results   : chr "asis"
##  $ echo      : logi FALSE
##  $ fig.height: num 5
##  $ fig.width : num 6
##  $ out.width : chr "80%"
##  $ comment   : logi NA
##  $ message   : logi FALSE
## 
## 
  |                                                                            
  |..................................................................    |  94%
##   ordinary text without R code
## 
## 
  |                                                                            
  |....................................................................  |  97%
## label: session_info (with options) 
## List of 3
##  $ include: logi TRUE
##  $ echo   : logi TRUE
##  $ results: chr "markup"
## 
## 
  |                                                                            
  |......................................................................| 100%
##   ordinary text without R code
## output file: /media/mdz/bfx6/bfx/adam_trewin/aln/mitch.knit.md
## /usr/bin/pandoc +RTS -K512m -RTS /media/mdz/bfx6/bfx/adam_trewin/aln/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpwUhoHt/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --standalone --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --include-in-header /tmp/RtmpwUhoHt/rmarkdown-str157dc2eb76f80.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'
## 
## Output created: /tmp/RtmpwUhoHt/mitch_report.html
## [1] TRUE

Session info

sessionInfo()
## R version 4.2.0 (2022-04-22)
## 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/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## 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] pkgload_1.2.4               GGally_2.1.2               
##  [3] ggplot2_3.3.5               beeswarm_0.4.0             
##  [5] gtools_3.9.2                tibble_3.1.6               
##  [7] dplyr_1.0.8                 echarts4r_0.4.3            
##  [9] kableExtra_1.3.4            gplots_3.1.1               
## [11] mitch_1.5.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.0           
## [21] BiocGenerics_0.38.0         reshape2_1.4.4             
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-3       ellipsis_0.3.2         rprojroot_2.0.2       
##  [4] XVector_0.32.0         rstudioapi_0.13        bit64_4.0.5           
##  [7] AnnotationDbi_1.54.1   fansi_1.0.2            xml2_1.3.3            
## [10] splines_4.2.0          cachem_1.0.6           geneplotter_1.70.0    
## [13] knitr_1.37             jsonlite_1.7.3         annotate_1.70.0       
## [16] png_0.1-7              shiny_1.7.1            compiler_4.2.0        
## [19] httr_1.4.2             assertthat_0.2.1       Matrix_1.4-1          
## [22] fastmap_1.1.0          cli_3.2.0              later_1.3.0           
## [25] htmltools_0.5.2        tools_4.2.0            gtable_0.3.0          
## [28] glue_1.6.1             GenomeInfoDbData_1.2.6 Rcpp_1.0.8            
## [31] jquerylib_0.1.4        vctrs_0.3.8            Biostrings_2.60.2     
## [34] svglite_2.1.0          xfun_0.29              stringr_1.4.0         
## [37] brio_1.1.3             testthat_3.1.2         rvest_1.0.2           
## [40] mime_0.12              lifecycle_1.0.1        XML_3.99-0.8          
## [43] zlibbioc_1.38.0        MASS_7.3-57            scales_1.1.1          
## [46] promises_1.2.0.1       RColorBrewer_1.1-2     yaml_2.3.5            
## [49] memoise_2.0.1          gridExtra_2.3          sass_0.4.0            
## [52] reshape_0.8.8          stringi_1.7.6          RSQLite_2.2.10        
## [55] highr_0.9              genefilter_1.74.0      desc_1.4.0            
## [58] caTools_1.18.2         BiocParallel_1.26.2    rlang_1.0.1           
## [61] pkgconfig_2.0.3        systemfonts_1.0.4      bitops_1.0-7          
## [64] evaluate_0.15          lattice_0.20-45        purrr_0.3.4           
## [67] htmlwidgets_1.5.4      bit_4.0.4              tidyselect_1.1.2      
## [70] plyr_1.8.6             magrittr_2.0.2         R6_2.5.1              
## [73] generics_0.1.2         DelayedArray_0.18.0    DBI_1.1.2             
## [76] pillar_1.7.0           withr_2.4.3            survival_3.2-13       
## [79] KEGGREST_1.32.0        RCurl_1.98-1.6         crayon_1.5.0          
## [82] KernSmooth_2.23-20     utf8_1.2.2             rmarkdown_2.11        
## [85] locfit_1.5-9.4         grid_4.2.0             blob_1.2.2            
## [88] digest_0.6.29          webshot_0.5.2          xtable_1.8-4          
## [91] httpuv_1.6.5           munsell_0.5.0          viridisLite_0.4.0     
## [94] bslib_0.3.1