Transcriptome sequencing of whole cell and mito fractions of sedentary and trained rats

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 (Ensembl version 99) with Kallisto then imported to R for analysis with DESeq2. Pathway level analysis was then done using mitch with Reactome gene sets.

Read counts

Import the Kallisto 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
## ENSRNOT00000000008.4   0.0000   5.17195  10.0000   6.12831   5.26551   3.0000
## ENSRNOT00000000009.5  10.0209  47.86210  53.7012  35.32830  34.89250  51.1897
## ENSRNOT00000000010.5   1.0000   1.00000   0.0000   0.00000   0.00000   0.0000
## ENSRNOT00000000011.5   0.0000   0.00000   1.0000   1.00000   0.00000   1.0000
## ENSRNOT00000000013.5 299.6020 674.72100 346.8700 553.52400 741.68800 466.5910
## ENSRNOT00000000018.7   0.0000   0.00000   2.0000   2.00000   1.00000   0.0000
##                             15       16      17        18       19        2
## ENSRNOT00000000008.4   3.90088   0.0000   0.000   4.07072  16.4054   1.3407
## ENSRNOT00000000009.5  49.16160  27.3429  83.000  51.24440  63.8893  23.2207
## ENSRNOT00000000010.5   0.00000   0.0000   0.000   0.00000   0.0000   0.0000
## ENSRNOT00000000011.5   0.00000   1.0000   1.000   0.00000   1.0000   4.0000
## ENSRNOT00000000013.5 341.29800 372.3020 416.936 509.46900 513.9120 320.8710
## ENSRNOT00000000018.7   0.00000   0.0000   3.000   0.00000   1.0000   1.0000
##                            20       21       22        23        24       25
## ENSRNOT00000000008.4   0.0000  20.9863   0.0000   6.54813   3.93027  13.1300
## ENSRNOT00000000009.5  42.3625  68.9222  17.6006  36.76790  24.36560  47.0854
## ENSRNOT00000000010.5   0.0000   2.0000   0.0000   0.00000   0.00000   0.0000
## ENSRNOT00000000011.5   2.0000   1.0000   6.0000   3.00000   3.00000   0.0000
## ENSRNOT00000000013.5 346.0990 936.4710 645.0050 845.56600 657.88000 653.9830
## ENSRNOT00000000018.7   0.0000   0.0000   0.0000   2.00000   4.00000   1.0000
##                             26       27       28       29        3        30
## ENSRNOT00000000008.4   3.92594  10.0000   0.0000   0.0000   7.0000   6.53318
## ENSRNOT00000000009.5  63.18030  76.1159  43.9129  42.8444  37.8331  54.25200
## ENSRNOT00000000010.5   0.00000   0.0000   0.0000   0.0000   1.0000   0.00000
## ENSRNOT00000000011.5   0.00000   0.0000   2.0000   2.0000   1.0000   0.00000
## ENSRNOT00000000013.5 699.78700 850.1850 771.5110 475.2300 491.6540 343.60600
## ENSRNOT00000000018.7   0.00000   0.0000   0.0000   0.0000   3.0000   0.00000
##                            31       32       33      34       35      36
## ENSRNOT00000000008.4   8.0000   0.0000  40.0898 19.1593  25.9218 11.4780
## ENSRNOT00000000009.5  41.2632  33.1978  28.8229 35.7135  70.6359 35.9501
## ENSRNOT00000000010.5   0.0000   1.0000  19.0000 11.0000  37.0000 15.0000
## ENSRNOT00000000011.5   4.0000   0.0000  35.0000 20.0000  34.0000 32.0000
## ENSRNOT00000000013.5 479.0410 319.9350 111.7250 72.2435 153.0400 62.8931
## ENSRNOT00000000018.7   0.0000   0.0000  14.0000 15.0000  29.0000 18.0000
##                           37      38      39         4       40      41
## ENSRNOT00000000008.4 13.7003 22.6135 14.9332   6.50394  8.52695 12.2591
## ENSRNOT00000000009.5 50.6873 54.7243 35.2301  27.46480 34.80920 26.7389
## ENSRNOT00000000010.5 40.0000 41.0000 35.0000   0.00000 29.00000 10.0000
## ENSRNOT00000000011.5 47.0000 33.0000 14.0000   0.00000 28.00000 22.0000
## ENSRNOT00000000013.5 52.7508 85.5584 68.9108 294.70900 55.07810 72.7735
## ENSRNOT00000000018.7 27.0000 17.0000 13.0000   0.00000 14.00000 15.0000
##                            42       43       44       45       46       47
## ENSRNOT00000000008.4  4.69899  13.6271  16.2942  24.1949  12.9460  32.2802
## ENSRNOT00000000009.5 33.73670 114.9440  47.4423  53.6637  41.5453  67.1163
## ENSRNOT00000000010.5 16.00000  35.0000  13.0000  32.0000  18.0000  25.0000
## ENSRNOT00000000011.5 21.00000  49.0000  29.0000  44.0000  39.0000  37.0000
## ENSRNOT00000000013.5 65.41440 121.1980 114.3950 131.3720 124.3370 131.1560
## ENSRNOT00000000018.7 11.00000  24.0000  15.0000  21.0000  18.0000  10.0000
##                            48      49        5       50       51      52
## ENSRNOT00000000008.4  16.3970 36.6810   0.0000  32.6092  20.9075 23.8737
## ENSRNOT00000000009.5  36.0636 49.3216  32.9113  53.3994  60.2464 36.7990
## ENSRNOT00000000010.5  20.0000 38.0000   0.0000  42.0000  46.0000 40.0000
## ENSRNOT00000000011.5  25.0000 43.0000   0.0000  48.0000  56.0000 33.0000
## ENSRNOT00000000013.5 114.0700 78.2890 824.7450 130.2330 150.1700 93.1834
## ENSRNOT00000000018.7  14.0000 25.0000   0.0000  22.0000  29.0000 17.0000
##                            53       54       55       56        57       58
## ENSRNOT00000000008.4  19.3599  24.1487  14.2184  32.6797   9.35056  10.1055
## ENSRNOT00000000009.5  39.9345  36.7259  56.0723  34.3157  55.51630  40.9847
## ENSRNOT00000000010.5  30.0000  40.0000  41.0000  29.0000  43.00000  27.0000
## ENSRNOT00000000011.5  54.0000  33.0000  53.0000  42.0000  33.00000  34.0000
## ENSRNOT00000000013.5 120.6350 129.7970 127.9770 122.6540 145.20200 117.3270
## ENSRNOT00000000018.7  18.0000  15.0000  24.0000  21.0000  17.00000  15.0000
##                           59        6       60       61       62       63
## ENSRNOT00000000008.4 31.4681  15.3007  17.6768  18.0824  21.5445  30.3799
## ENSRNOT00000000009.5 41.2987  34.8236  53.4337  70.0980  50.4952  49.1870
## ENSRNOT00000000010.5 24.0000   0.0000  24.0000  34.0000  33.0000  41.0000
## ENSRNOT00000000011.5 34.0000   0.0000  50.0000  47.0000  45.0000  48.0000
## ENSRNOT00000000013.5 96.5047 645.2290 100.4960 139.6330 153.3700 129.0280
## ENSRNOT00000000018.7 14.0000   0.0000  25.0000  14.0000  19.0000  25.0000
##                            64         7        8        9 95      96
## ENSRNOT00000000008.4  29.0847   1.67668   0.0000   0.0000  0  0.0000
## ENSRNOT00000000009.5  48.8472  19.33470  22.7956  30.2034  0  0.0000
## ENSRNOT00000000010.5  24.0000   0.00000   1.0000   0.0000  0  0.0000
## ENSRNOT00000000011.5  51.0000   0.00000   3.0000   0.0000  0  0.0000
## ENSRNOT00000000013.5 155.7100 400.78800 316.4870 461.2840  6 28.8082
## ENSRNOT00000000018.7  33.0000   0.00000   3.0000   0.0000  0  0.0000

#dont forget gene names
g<-read.table("../ref/Rattus_norvegicus.Rnor_6.0.cdna+ncrna.gene_names.tsv",row.names=1)
g$gene_ID <- paste(g$V2,g$V3,g$V4,g$V5)
head(g)

##                      V2                   V3      V4             V5
## ENSRNOT00000047550.4 MT ENSRNOG00000030644.4  Mt-nd1 protein_coding
## ENSRNOT00000040993.4 MT ENSRNOG00000031033.4  Mt-nd2 protein_coding
## ENSRNOT00000050156.3 MT ENSRNOG00000034234.3  Mt-co1 protein_coding
## ENSRNOT00000043693.3 MT ENSRNOG00000030371.3  Mt-co2 protein_coding
## ENSRNOT00000046201.3 MT ENSRNOG00000033299.3 Mt-atp8 protein_coding
## ENSRNOT00000046108.3 MT ENSRNOG00000031979.3 Mt-atp6 protein_coding
##                                                             gene_ID
## ENSRNOT00000047550.4  MT ENSRNOG00000030644.4 Mt-nd1 protein_coding
## ENSRNOT00000040993.4  MT ENSRNOG00000031033.4 Mt-nd2 protein_coding
## ENSRNOT00000050156.3  MT ENSRNOG00000034234.3 Mt-co1 protein_coding
## ENSRNOT00000043693.3  MT ENSRNOG00000030371.3 Mt-co2 protein_coding
## ENSRNOT00000046201.3 MT ENSRNOG00000033299.3 Mt-atp8 protein_coding
## ENSRNOT00000046108.3 MT ENSRNOG00000031979.3 Mt-atp6 protein_coding

g <- g[,ncol(g),drop=FALSE]
x<-merge(g,x,by=0)
rownames(x) <- x[,1]
x[,1]=NULL
# aggregate Tx data to genes
xx <- aggregate(. ~ gene_ID,x,sum)
# now round to integers so that DESeq2 doesn't fail
rownames(xx) <- xx[,1]
xx[,1]=NULL
x <- round(xx)
head(x)

##                                                       1   10   11   12   13
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        3301 4795 4460 4965 5104
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding    0    0    0    0    0
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding        1525 1671 1845 1952 1631
## 1 ENSRNOG00000001489.5 Tbp protein_coding           189  181  195  156  154
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding          95  114  123  125  120
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding        188  150  155  192  199
##                                                      14   15   16   17   18
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        4177 3716 4657 4314 4337
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding    0    0    0    1    0
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding        1286 1407 1245 1756 1851
## 1 ENSRNOG00000001489.5 Tbp protein_coding           156  120  159  188  197
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding          93  107   96  190  159
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding        173  158  192  250  262
##                                                      19    2   20   21   22
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        3901 2782 3403 5194 3945
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding    0    0    0    0    0
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding        1947 1204 2112 2717 1679
## 1 ENSRNOG00000001489.5 Tbp protein_coding           210  147  165  241  183
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding         175   59  171  227  124
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding        188  106  210  365  253
##                                                      23   24   25   26   27
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        4823 4317 5347 4593 5806
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding    0    0    0    0    0
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding        1861 1707 1830 1644 1567
## 1 ENSRNOG00000001489.5 Tbp protein_coding           234  200  201  194  153
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding         137  152  177  118  129
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding        246  323  300  266  194
##                                                      28   29    3   30   31
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        4751 3516 3178 3746 3497
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding    0    0    0    0    0
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding        1776 1343 1737 1470 1762
## 1 ENSRNOG00000001489.5 Tbp protein_coding           155  135  127  160  171
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding         161   62  143  112  123
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding        248  123  110  217  155
##                                                      32  33  34  35  36  37  38
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        2733 192 137 252 173 206 219
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding    0  20  24  23  13  15  19
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding        1083  27  15  30  28  16  29
## 1 ENSRNOG00000001489.5 Tbp protein_coding           128  19  23  31  18  28  45
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding          84  21  12  33  12  13  30
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding        183 250 238 360 245 335 349
##                                                     39    4  40  41  42  43  44
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        166 3357 168 154 144 336 190
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding  24    0  15  18  13  17   8
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding         23 1649  23  18  11  43  19
## 1 ENSRNOG00000001489.5 Tbp protein_coding           39  168  17  25  21  54  26
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding         20   85  19  13  14  31  28
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding       229  185 227 241 191 482 286
##                                                     45  46  47  48  49    5  50
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        188 199 205 162 229 4574 243
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding  25  12  10  17  21    0  31
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding         23  28  22  21  20 1879  37
## 1 ENSRNOG00000001489.5 Tbp protein_coding           38  25  34  24  30  150  43
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding         18  24  18  25  19  128  29
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding       323 305 339 250 340  151 378
##                                                     51  52  53  54  55  56  57
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        296 243 214 232 249 265 237
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding  44  11  26  28  29  37  32
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding         40  24  24  29  29  22  29
## 1 ENSRNOG00000001489.5 Tbp protein_coding           42  36  44  45  41  44  26
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding         28  35  26  33  28  30  20
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding       501 384 364 336 355 420 408
##                                                     58  59    6  60  61  62  63
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        235 185 4385 281 255 250 300
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding  14  30    0  24  24  31  31
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding         43  34 1459  45  25  37  34
## 1 ENSRNOG00000001489.5 Tbp protein_coding           36  25  150  37  34  46  42
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding         20  19  105  24  35  25  28
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding       403 345  183 382 442 448 440
##                                                     64    7    8    9 95 96
## 1 ENSRNOG00000000417.7 Numa1 protein_coding        247 3642 3364 3272  2 32
## 1 ENSRNOG00000001466.6 LOC100361492 protein_coding  23    0    1    0  0  0
## 1 ENSRNOG00000001488.6 Psmb1 protein_coding         19 1518 1201 1268  5  0
## 1 ENSRNOG00000001489.5 Tbp protein_coding           25  178   88  143  0  0
## 1 ENSRNOG00000001490.4 Pdcd2 protein_coding         22  108   66   78  0  0
## 1 ENSRNOG00000001492.7 Slc8a2 protein_coding       419  147  164  127  0  0

write.table(x,file="countmatrix_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))
mtfrac <-  colSums(x[grep("^MT",rownames(x)),]) / 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.3338761

median(wholeT)

## [1] 0.3550418

table(mitoS>0.4)

## 
## FALSE  TRUE 
##    14     2

table(mitoT>0.4)

## 
## FALSE  TRUE 
##    11     5

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	21258507	0.2601233
10	10	wholetissue	T	10	wholetissue T	32402582	0.3723191
11	11	wholetissue	S	11	wholetissue S	27624509	0.3521941
12	12	wholetissue	S	12	wholetissue S	29209069	0.2845851
13	13	wholetissue	T	13	wholetissue T	29174764	0.3195613
14	14	wholetissue	T	14	wholetissue T	26930441	0.3552110
15	15	wholetissue	S	15	wholetissue S	24139589	0.3354321
16	16	wholetissue	S	16	wholetissue S	26585363	0.3323202
17	17	wholetissue	T	17	wholetissue T	29764939	0.3617978
18	18	wholetissue	T	18	wholetissue T	32286871	0.3904784
19	19	wholetissue	S	19	wholetissue S	29474137	0.3574090
2	2	wholetissue	T	2	wholetissue T	18839052	0.3441832
20	20	wholetissue	S	20	wholetissue S	29473669	0.3501682
21	21	wholetissue	T	21	wholetissue T	38356697	0.3680570
22	22	wholetissue	T	22	wholetissue T	26332319	0.3551245
23	23	wholetissue	S	23	wholetissue S	31385455	0.3802950
24	24	wholetissue	S	24	wholetissue S	27209791	0.2678982
25	25	wholetissue	T	25	wholetissue T	33288494	0.3205449
26	26	wholetissue	T	26	wholetissue T	28504502	0.2691352
27	27	wholetissue	T	27	wholetissue T	30196627	0.3549592
28	28	wholetissue	T	28	wholetissue T	31275863	0.2879719
29	29	wholetissue	S	29	wholetissue S	24900845	0.3670548
3	3	wholetissue	S	3	wholetissue S	23495593	0.3273260
30	30	wholetissue	S	30	wholetissue S	25393464	0.3000018
31	31	wholetissue	S	31	wholetissue S	26748801	0.3738198
32	32	wholetissue	S	32	wholetissue S	19491222	0.3556323
33	33	mito	T	1	mito T	7938936	0.6465534
34	34	mito	T	2	mito T	12294342	0.8147126
35	35	mito	S	3	mito S	3781587	0.1980579
36	36	mito	S	4	mito S	2752655	0.2414353
37	37	mito	T	5	mito T	3445116	0.2070290
38	38	mito	T	6	mito T	3101202	0.0577676
39	39	mito	S	7	mito S	3516550	0.3757165
4	4	wholetissue	S	4	wholetissue S	23993574	0.2853209
40	40	mito	S	8	mito S	2164670	0.1035040
41	41	mito	T	9	mito T	2041406	0.0597436
42	42	mito	T	10	mito T	1731503	0.0656817
43	43	mito	S	11	mito S	4641993	0.1646653
44	44	mito	S	12	mito S	2576114	0.0634708
45	45	mito	T	13	mito T	3027547	0.0948844
46	46	mito	T	14	mito T	3000299	0.1706080
47	47	mito	S	15	mito S	11142807	0.7268125
48	48	mito	S	16	mito S	3010541	0.2668693
49	49	mito	T	17	mito T	3080929	0.1054198
5	5	wholetissue	T	5	wholetissue T	30674155	0.3747933
50	50	mito	T	18	mito T	4044280	0.1292304
51	51	mito	S	19	mito S	4890168	0.1009622
52	52	mito	S	20	mito S	6589045	0.3775514
53	53	mito	T	21	mito T	4720580	0.2118263
54	54	mito	T	22	mito T	3685725	0.1458188
55	55	mito	S	23	mito S	3101674	0.0446639
56	56	mito	S	24	mito S	5790497	0.2490106
57	57	mito	T	25	mito T	4620986	0.2491953
58	58	mito	T	26	mito T	6786437	0.4821588
59	59	mito	T	27	mito T	10802671	0.6646957
6	6	wholetissue	T	6	wholetissue T	25369115	0.3510836
60	60	mito	T	28	mito T	9524968	0.5887804
61	61	mito	S	29	mito S	4721621	0.2121115
62	62	mito	S	30	mito S	4409979	0.1527370
63	63	mito	S	31	mito S	4277756	0.1687420
64	64	mito	S	32	mito S	7138202	0.5016151
7	7	wholetissue	S	7	wholetissue S	24756089	0.3214936
8	8	wholetissue	S	8	wholetissue S	20318569	0.2543348
9	9	wholetissue	T	9	wholetissue T	24443981	0.3754423

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")

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")

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 207 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

##  chr [1:21] "9 ENSRNOG00000019961.5 Tcte1 protein_coding" ...
##  chr [1:2] "19 ENSRNOG00000015519.6 Ces1d protein_coding" ...

as.data.frame(de1[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
9 ENSRNOG00000019961.5 Tcte1 protein_coding	63.23168	0.9887063	0.1472459	6.714661	0.00e+00	0.0000003
9 ENSRNOG00000019958.5 Tmem151b protein_coding	43.76548	1.5367206	0.2476454	6.205326	0.00e+00	0.0000044
4 ENSRNOG00000060514.1 LOC103692166 protein_coding	41.29021	0.7035412	0.1333954	5.274104	1.00e-07	0.0007247
8 ENSRNOG00000015896.7 Rbpms2 protein_coding	64.48790	0.5181886	0.1077025	4.811297	1.50e-06	0.0061091
3 ENSRNOG00000008482.7 Rbms1 protein_coding	485.10997	0.1890229	0.0399331	4.733485	2.20e-06	0.0071930
19 ENSRNOG00000012482.7 Ndrg4 protein_coding	119.21510	0.7977010	0.1698899	4.695400	2.70e-06	0.0072269
14 ENSRNOG00000000060.7 Cplx1 protein_coding	18.93501	0.9098606	0.1965642	4.628822	3.70e-06	0.0085613
4 ENSRNOG00000010284.6 St3gal5 protein_coding	416.51249	0.4223257	0.0918658	4.597202	4.30e-06	0.0087225
1 ENSRNOG00000015071.4 Zim1 protein_coding	520.46906	0.4077851	0.0933999	4.366012	1.27e-05	0.0221628
13 ENSRNOG00000007290.6 Atp1a2 protein_coding	37352.46595	0.1995007	0.0458599	4.350219	1.36e-05	0.0221628
3 ENSRNOG00000014710.6 Prom2 protein_coding	154.64455	0.4353713	0.1011122	4.305824	1.66e-05	0.0246463
13 ENSRNOG00000005093.4 Lgr6 protein_coding	74.37305	0.6745253	0.1593156	4.233895	2.30e-05	0.0288562
19 ENSRNOG00000015519.6 Ces1d protein_coding	163.95783	-0.6050191	0.1429161	-4.233387	2.30e-05	0.0288562
7 ENSRNOG00000006086.7 Lynx1 protein_coding	1233.69841	0.2335705	0.0561119	4.162585	3.15e-05	0.0357858
12 ENSRNOG00000001115.5 Slc29a4 protein_coding	16.27318	1.1606653	0.2795351	4.152128	3.29e-05	0.0357858
5 ENSRNOG00000023812.5 Raver2 protein_coding	112.30813	0.4713947	0.1141614	4.129196	3.64e-05	0.0363784
10 ENSRNOG00000046195.2 AABR07029596.1 protein_coding	50.30723	0.7663390	0.1861465	4.116860	3.84e-05	0.0363784
15 ENSRNOG00000011323.5 Lgi3 protein_coding	157.83997	0.4762472	0.1159760	4.106429	4.02e-05	0.0363784
10 ENSRNOG00000009370.5 Tbkbp1 protein_coding	232.31212	0.3928446	0.0959925	4.092452	4.27e-05	0.0366090
1 ENSRNOG00000016356.7 Got1 protein_coding	9196.33671	0.4218979	0.1034650	4.077687	4.55e-05	0.0370619

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 43 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
11 ENSRNOG00000057596.1 AABR07033056.1 processed_pseudogene	6.799236	-0.9858589	0.2378721	-4.144492	0.0000341	0.8164156
20 ENSRNOG00000033652.3 AABR07045321.1 protein_coding	31.385811	0.4180776	0.1107243	3.775844	0.0001595	0.8164156
16 ENSRNOG00000031335.4 Ankrd37 protein_coding	20.342893	0.5214769	0.1381419	3.774936	0.0001600	0.8164156
17 ENSRNOG00000018698.7 Wac protein_coding	69.447049	0.2927467	0.0784855	3.729946	0.0001915	0.8164156
7 ENSRNOG00000053181.1 Aqp6 protein_coding	28.839273	-0.4517375	0.1211477	-3.728817	0.0001924	0.8164156
10 ENSRNOG00000002388.4 Trim41 protein_coding	72.206861	-0.2721179	0.0744236	-3.656338	0.0002558	0.8164156
2 ENSRNOG00000058603.1 AABR07011434.1 lincRNA	83.746925	0.2713824	0.0748666	3.624881	0.0002891	0.8164156
8 ENSRNOG00000010196.5 Glb1 protein_coding	105.047842	-0.2683572	0.0742367	-3.614887	0.0003005	0.8164156
2 ENSRNOG00000059604.1 Lor protein_coding	37.381077	0.4132954	0.1145260	3.608747	0.0003077	0.8164156
12 ENSRNOG00000057927.1 AABR07035796.1 lincRNA	64.919074	-0.3337369	0.0927046	-3.600003	0.0003182	0.8164156
1 ENSRNOG00000013426.7 Mrgprf protein_coding	68.012764	-0.3029277	0.0843878	-3.589710	0.0003310	0.8164156
18 ENSRNOG00000057832.1 Rnf125 protein_coding	24.503842	0.5066023	0.1414879	3.580535	0.0003429	0.8164156
4 ENSRNOG00000000021.3 AABR07061902.1 pseudogene	31.184421	-0.4398168	0.1230320	-3.574815	0.0003505	0.8164156
2 ENSRNOG00000051623.1 Exosc8 protein_coding	22.769240	0.4874604	0.1374934	3.545336	0.0003921	0.8164156
1 ENSRNOG00000043320.3 LOC103691092 protein_coding	36.149436	0.4002888	0.1133295	3.532080	0.0004123	0.8164156
19 ENSRNOG00000011635.7 Ces2e protein_coding	50.953991	0.3319378	0.0940948	3.527697	0.0004192	0.8164156
14 ENSRNOG00000002271.6 Slain2 protein_coding	70.124559	0.2649451	0.0752231	3.522125	0.0004281	0.8164156
10 ENSRNOG00000002412.7 Mapk7 protein_coding	54.211409	-0.4041696	0.1169384	-3.456261	0.0005477	0.8503168
1 ENSRNOG00000057034.1 AABR07003841.1 lincRNA	9.142742	1.0409635	0.3019970	3.446934	0.0005670	0.8503168
13 ENSRNOG00000022724.5 Cnih3 protein_coding	18.155512	-0.4739114	0.1379474	-3.435449	0.0005916	0.8503168

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 205 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

##  chr [1:22549] "1 ENSRNOG00000010983.6 Otog protein_coding" ...
##  chr [1:7732] "1 ENSRNOG00000000417.7 Numa1 protein_coding" ...

as.data.frame(de3[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
1 ENSRNOG00000000417.7 Numa1 protein_coding	1224.2153	-2.278258	0.0515175	-44.22298	0	0
1 ENSRNOG00000001488.6 Psmb1 protein_coding	451.8783	-4.128107	0.0938347	-43.99342	0	0
1 ENSRNOG00000001494.6 Napa protein_coding	445.0417	-3.135309	0.0695501	-45.07986	0	0
1 ENSRNOG00000010117.7 Eif3a protein_coding	1420.1411	-4.196810	0.0704317	-59.58695	0	0
1 ENSRNOG00000010427.5 Ipo7 protein_coding	628.0483	-2.832880	0.0647304	-43.76431	0	0
1 ENSRNOG00000010958.5 Prdx3 protein_coding	609.2282	-4.933490	0.1091842	-45.18500	0	0
1 ENSRNOG00000010964.6 Akap13 protein_coding	1337.5382	-2.429310	0.0552866	-43.94034	0	0
1 ENSRNOG00000010983.6 Otog protein_coding	429.4795	6.229600	0.1128381	55.20832	0	0
1 ENSRNOG00000011282.6 Ppp2r1a protein_coding	547.8787	-3.099354	0.0772535	-40.11927	0	0
1 ENSRNOG00000011339.6 Slk protein_coding	411.9456	-2.708499	0.0651728	-41.55871	0	0
1 ENSRNOG00000011346.5 Ehd2 protein_coding	787.3396	-4.212573	0.1028450	-40.96042	0	0
1 ENSRNOG00000011460.6 Arfgef3 protein_coding	126.4895	6.147113	0.1524340	40.32640	0	0
1 ENSRNOG00000011659.5 Alpk3 protein_coding	4290.1696	-5.087543	0.0763225	-66.65851	0	0
1 ENSRNOG00000011885.5 Rhpn2 protein_coding	311.6309	3.855666	0.0936647	41.16456	0	0
1 ENSRNOG00000011945.7 Cyfip1 protein_coding	412.5033	-1.842031	0.0458336	-40.18952	0	0
1 ENSRNOG00000012110.7 Col17a1 protein_coding	162.0319	5.013767	0.0999438	50.16584	0	0
1 ENSRNOG00000012142.6 Glyat protein_coding	297.2643	4.700972	0.1054415	44.58370	0	0
1 ENSRNOG00000012206.6 Dennd5a protein_coding	601.6843	-2.443968	0.0460011	-53.12842	0	0
1 ENSRNOG00000012383.6 Ndufc2 protein_coding	655.0776	-5.849854	0.1420982	-41.16770	0	0
1 ENSRNOG00000012415.6 Mpc1 protein_coding	1625.1625	-6.426814	0.1195915	-53.73971	0	0

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 105 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

##  chr [1:22544] "1 ENSRNOG00000009932.8 Mctp2 protein_coding" ...
##  chr [1:7674] "1 ENSRNOG00000000417.7 Numa1 protein_coding" ...

as.data.frame(de4[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
1 ENSRNOG00000000417.7 Numa1 protein_coding	1262.9612	-2.385122	0.0462115	-51.61316	0	0
1 ENSRNOG00000001488.6 Psmb1 protein_coding	427.6390	-3.938947	0.0933781	-42.18278	0	0
1 ENSRNOG00000001494.6 Napa protein_coding	427.3646	-3.115259	0.0788663	-39.50050	0	0
1 ENSRNOG00000009932.8 Mctp2 protein_coding	138.4064	3.880420	0.0862188	45.00668	0	0
1 ENSRNOG00000010117.7 Eif3a protein_coding	1418.0006	-4.172418	0.0872323	-47.83111	0	0
1 ENSRNOG00000010427.5 Ipo7 protein_coding	627.3105	-2.669449	0.0709960	-37.60002	0	0
1 ENSRNOG00000010958.5 Prdx3 protein_coding	621.0700	-4.616300	0.1033410	-44.67054	0	0
1 ENSRNOG00000010964.6 Akap13 protein_coding	1329.1460	-2.320667	0.0496862	-46.70648	0	0
1 ENSRNOG00000010983.6 Otog protein_coding	401.5225	6.040622	0.1117856	54.03755	0	0
1 ENSRNOG00000011058.7 Utrn protein_coding	1084.9255	-1.754523	0.0463132	-37.88383	0	0
1 ENSRNOG00000011339.6 Slk protein_coding	421.8244	-2.587643	0.0610473	-42.38754	0	0
1 ENSRNOG00000011346.5 Ehd2 protein_coding	897.6189	-4.372625	0.0938367	-46.59823	0	0
1 ENSRNOG00000011460.6 Arfgef3 protein_coding	118.8150	6.636153	0.1765664	37.58445	0	0
1 ENSRNOG00000011659.5 Alpk3 protein_coding	4371.2836	-5.090869	0.0738495	-68.93577	0	0
1 ENSRNOG00000011885.5 Rhpn2 protein_coding	311.8998	3.750411	0.0966640	38.79844	0	0
1 ENSRNOG00000012110.7 Col17a1 protein_coding	158.8087	5.136959	0.1051545	48.85152	0	0
1 ENSRNOG00000012142.6 Glyat protein_coding	298.2102	4.774947	0.1086791	43.93622	0	0
1 ENSRNOG00000012206.6 Dennd5a protein_coding	623.9719	-2.682520	0.0518118	-51.77429	0	0
1 ENSRNOG00000012383.6 Ndufc2 protein_coding	635.1984	-5.599595	0.1311112	-42.70875	0	0
1 ENSRNOG00000012415.6 Mpc1 protein_coding	1629.2379	-6.332387	0.1117885	-56.64613	0	0

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.

2. Mito fraction S versus T: 0 DEGs

3. Whole versus mito fraction in S: 20 DEGs

4. Whole versus mito fraction in T: 0 DEGs

5. 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	7938936	0.6465534	1
34	34	mito	T	2	mito T	12294342	0.8147126	1
47	47	mito	S	15	mito S	11142807	0.7268125	0
58	58	mito	T	26	mito T	6786437	0.4821588	1
59	59	mito	T	27	mito T	10802671	0.6646957	1
60	60	mito	T	28	mito T	9524968	0.5887804	1
64	64	mito	S	32	mito S	7138202	0.5016151	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
17 ENSRNOG00000056825.2 AABR07072539.4 processed_pseudogene	11.661670	6.5676953	1.5385473	4.268764	0.0000197	0.6313871
1 ENSRNOG00000020025.4 LOC108348052 protein_coding	17.782509	2.9172245	0.7530693	3.873780	0.0001072	0.9999479
3 ENSRNOG00000050864.2 Cpne1 protein_coding	14.072847	-2.9272763	0.7602602	-3.850361	0.0001179	0.9999479
2 ENSRNOG00000020832.4 C2cd4d protein_coding	15.822890	2.7222352	0.7431078	3.663311	0.0002490	0.9999479
7 ENSRNOG00000050884.2 AABR07058745.1 protein_coding	50.195404	1.6687089	0.4779421	3.491445	0.0004804	0.9999479
20 ENSRNOG00000045654.2 LOC108348108 protein_coding	7.187484	5.8570854	1.6931454	3.459293	0.0005416	0.9999479
4 ENSRNOG00000049639.2 LOC100909784 protein_coding	52.027683	1.1890107	0.3485934	3.410881	0.0006475	0.9999479
13 ENSRNOG00000056279.1 AC115369.1 lincRNA	199.904660	-0.5561194	0.1700077	-3.271142	0.0010711	0.9999479
6 ENSRNOG00000051280.1 AABR07064502.1 processed_pseudogene	5.509108	5.4622037	1.6700748	3.270634	0.0010731	0.9999479
12 ENSRNOG00000057927.1 AABR07035796.1 lincRNA	58.135739	-0.9607682	0.2975668	-3.228748	0.0012433	0.9999479
18 ENSRNOG00000051569.1 LOC103690147 protein_coding	16.613644	2.4049396	0.7510958	3.201908	0.0013652	0.9999479
14 ENSRNOG00000050406.1 AABR07014914.1 miRNA	102.689535	-0.7408943	0.2335580	-3.172207	0.0015129	0.9999479
5 ENSRNOG00000048618.1 AABR07049890.1 miRNA	68.590808	0.9768229	0.3118632	3.132216	0.0017349	0.9999479
2 ENSRNOG00000061049.1 AABR07011512.1 lincRNA	391.297383	0.4903412	0.1577683	3.107983	0.0018837	0.9999479
1 ENSRNOG00000038697.3 AABR07003616.1 protein_coding	2.349799	-4.7978937	1.5525918	-3.090248	0.0019999	0.9999479
MT ENSRNOG00000029070.3 AY172581.1 Mt_tRNA	28.057647	1.9396338	0.6302285	3.077667	0.0020863	0.9999479
17 ENSRNOG00000052667.1 AABR07027555.1 processed_pseudogene	4.412753	5.1807393	1.6974583	3.052057	0.0022728	0.9999479
KL568337.1 ENSRNOG00000057871.1 5_8S_rRNA rRNA	5.377655	-3.1465703	1.0320391	-3.048887	0.0022969	0.9999479
13 ENSRNOG00000052639.1 LOC108352591 lincRNA	11.174558	2.3679557	0.7780118	3.043599	0.0023377	0.9999479
8 ENSRNOG00000059264.1 U6 snRNA	5.160780	5.3206380	1.7779022	2.992649	0.0027657	0.9999479

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	27624509	0.3521941	0
12	12	wholetissue	S	12	wholetissue S	29209069	0.2845851	0
15	15	wholetissue	S	15	wholetissue S	24139589	0.3354321	0
16	16	wholetissue	S	16	wholetissue S	26585363	0.3323202	0
19	19	wholetissue	S	19	wholetissue S	29474137	0.3574090	0
20	20	wholetissue	S	20	wholetissue S	29473669	0.3501682	0
23	23	wholetissue	S	23	wholetissue S	31385455	0.3802950	0
24	24	wholetissue	S	24	wholetissue S	27209791	0.2678982	0
29	29	wholetissue	S	29	wholetissue S	24900845	0.3670548	0
3	3	wholetissue	S	3	wholetissue S	23495593	0.3273260	0
30	30	wholetissue	S	30	wholetissue S	25393464	0.3000018	0
31	31	wholetissue	S	31	wholetissue S	26748801	0.3738198	0
32	32	wholetissue	S	32	wholetissue S	19491222	0.3556323	0
4	4	wholetissue	S	4	wholetissue S	23993574	0.2853209	0
7	7	wholetissue	S	7	wholetissue S	24756089	0.3214936	0
8	8	wholetissue	S	8	wholetissue S	20318569	0.2543348	0
47	47	mito	S	15	mito S	11142807	0.7268125	1
64	64	mito	S	32	mito S	7138202	0.5016151	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 5629 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

##  chr [1:18] "9 ENSRNOG00000046834.3 C3 protein_coding" ...
##  chr [1:2] "14 ENSRNOG00000027791.5 AABR07015346.1 processed_pseudogene" ...

as.data.frame(de3[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
9 ENSRNOG00000046834.3 C3 protein_coding	189.947150	6.694343	0.8044297	8.321849	0.00e+00	0.0000000
14 ENSRNOG00000027791.5 AABR07015346.1 processed_pseudogene	154.480953	-23.408177	3.3688533	-6.948411	0.00e+00	0.0000000
13 ENSRNOG00000004405.5 Pigr protein_coding	33.147013	5.576134	0.8340084	6.685944	0.00e+00	0.0000002
2 ENSRNOG00000011200.3 Bhmt protein_coding	23.212978	4.704729	0.7534309	6.244406	0.00e+00	0.0000026
19 ENSRNOG00000014964.8 Hp protein_coding	85.923541	3.648683	0.6437731	5.667654	0.00e+00	0.0000625
16 ENSRNOG00000017381.8 Itih4 protein_coding	27.994702	6.468558	1.1437031	5.655802	0.00e+00	0.0000625
6 ENSRNOG00000005542.8 Apob protein_coding	91.699635	9.580368	1.7184669	5.574950	0.00e+00	0.0000840
19 ENSRNOG00000015438.5 LOC501233 protein_coding	19.393163	5.876676	1.0579672	5.554687	0.00e+00	0.0000840
4 ENSRNOG00000006709.6 Pzp protein_coding	51.516746	7.538147	1.3639175	5.526835	0.00e+00	0.0000875
1 ENSRNOG00000017223.7 Plg protein_coding	27.573652	6.470047	1.2003796	5.390001	1.00e-07	0.0001703
1 ENSRNOG00000046727.2 Abcc2 protein_coding	36.686475	3.034077	0.5814550	5.218078	2.00e-07	0.0003850
6 ENSRNOG00000010478.7 LOC299282 protein_coding	196.409234	2.161695	0.4150939	5.207726	2.00e-07	0.0003850
2 ENSRNOG00000012436.6 Adh6 protein_coding	12.016400	3.126573	0.6145946	5.087212	4.00e-07	0.0006754
1 ENSRNOG00000024016.7 Cyp2c6v1 protein_coding	19.504106	2.998256	0.6171350	4.858347	1.20e-06	0.0020431
3 ENSRNOG00000018899.8 C5 protein_coding	6.163963	5.651434	1.2054685	4.688164	2.80e-06	0.0044408
3 ENSRNOG00000010128.7 Slc27a2 protein_coding	6.107139	4.339468	0.9795889	4.429887	9.40e-06	0.0142390
4 ENSRNOG00000006675.5 Fabp1 protein_coding	52.467979	1.634886	0.3735441	4.376687	1.20e-05	0.0171276
1 ENSRNOG00000031391.4 Ceacam16 protein_coding	40.661418	-2.349823	0.5397098	-4.353863	1.34e-05	0.0179564
16 ENSRNOG00000033386.5 Itih1 protein_coding	11.068834	8.014672	1.8705981	4.284550	1.83e-05	0.0232876
2 ENSRNOG00000013736.6 C9 protein_coding	11.301353	3.221508	0.7727049	4.169131	3.06e-05	0.0369424

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	21258507	0.2601233	0
10	10	wholetissue	T	10	wholetissue T	32402582	0.3723191	0
13	13	wholetissue	T	13	wholetissue T	29174764	0.3195613	0
14	14	wholetissue	T	14	wholetissue T	26930441	0.3552110	0
17	17	wholetissue	T	17	wholetissue T	29764939	0.3617978	0
18	18	wholetissue	T	18	wholetissue T	32286871	0.3904784	0
2	2	wholetissue	T	2	wholetissue T	18839052	0.3441832	0
21	21	wholetissue	T	21	wholetissue T	38356697	0.3680570	0
22	22	wholetissue	T	22	wholetissue T	26332319	0.3551245	0
25	25	wholetissue	T	25	wholetissue T	33288494	0.3205449	0
26	26	wholetissue	T	26	wholetissue T	28504502	0.2691352	0
27	27	wholetissue	T	27	wholetissue T	30196627	0.3549592	0
28	28	wholetissue	T	28	wholetissue T	31275863	0.2879719	0
5	5	wholetissue	T	5	wholetissue T	30674155	0.3747933	0
6	6	wholetissue	T	6	wholetissue T	25369115	0.3510836	0
9	9	wholetissue	T	9	wholetissue T	24443981	0.3754423	0
33	33	mito	T	1	mito T	7938936	0.6465534	1
34	34	mito	T	2	mito T	12294342	0.8147126	1
58	58	mito	T	26	mito T	6786437	0.4821588	1
59	59	mito	T	27	mito T	10802671	0.6646957	1
60	60	mito	T	28	mito T	9524968	0.5887804	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 294 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(de4[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
17 ENSRNOG00000053995.1 AABR07027810.4 miRNA	8.350737	3.0222524	0.6855928	4.408232	1.04e-05	0.0900534
2 ENSRNOG00000048082.2 AABR07013776.1 protein_coding	8.905297	3.5112469	0.8304255	4.228250	2.36e-05	0.0900534
5 ENSRNOG00000055453.1 AC127935.1 snoRNA	3.827934	2.8729506	0.6804509	4.222128	2.42e-05	0.0900534
3 ENSRNOG00000011715.3 AABR07053669.1 pseudogene	9.923043	2.6658784	0.6323049	4.216128	2.49e-05	0.0900534
2 ENSRNOG00000045639.2 RGD1560883 protein_coding	70.610633	1.1875298	0.2836711	4.186291	2.84e-05	0.0900534
8 ENSRNOG00000023809.5 Opcml protein_coding	59.018462	1.0449219	0.2506728	4.168469	3.07e-05	0.0900534
1 ENSRNOG00000060312.1 LOC108348197 protein_coding	271.784529	-2.4832702	0.5972205	-4.158046	3.21e-05	0.0900534
1 ENSRNOG00000015717.8 Ptpre protein_coding	86.690130	0.9615012	0.2316621	4.150447	3.32e-05	0.0900534
7 ENSRNOG00000005551.6 Derl1 protein_coding	541.437086	0.2653631	0.0642131	4.132542	3.59e-05	0.0900534
5 ENSRNOG00000022953.4 Ccdc163 protein_coding	36.357350	1.6456259	0.3995540	4.118657	3.81e-05	0.0900534
3 ENSRNOG00000050251.1 MGC105649 protein_coding	13.324704	2.2234732	0.5398584	4.118623	3.81e-05	0.0900534
16 ENSRNOG00000045643.3 AABR07026271.2 protein_coding	5.871765	5.9098751	1.4396072	4.105200	4.04e-05	0.0900534
15 ENSRNOG00000047211.1 Fzd3 protein_coding	42.068994	1.2291144	0.3004629	4.090736	4.30e-05	0.0900534
9 ENSRNOG00000050660.1 LOC100911713 protein_coding	86.887371	-2.2825034	0.5595076	-4.079486	4.51e-05	0.0900534
13 ENSRNOG00000058096.1 AABR07021704.1 protein_coding	181.540160	0.6374446	0.1563429	4.077220	4.56e-05	0.0900534
X ENSRNOG00000002413.6 Gpc4 protein_coding	905.969267	1.2337233	0.3037032	4.062266	4.86e-05	0.0900534
4 ENSRNOG00000006726.7 Zfp9 protein_coding	63.475641	1.7474409	0.4332597	4.033241	5.50e-05	0.0928726
1 ENSRNOG00000043141.3 Ap3s2 protein_coding	109.562150	-0.6388746	0.1586301	-4.027450	5.64e-05	0.0928726
3 ENSRNOG00000016488.6 Pltp protein_coding	212.517369	-0.5918376	0.1475718	-4.010507	6.06e-05	0.0945435
2 ENSRNOG00000053898.1 AABR07013208.1 snoRNA	8.556559	1.7890899	0.4494923	3.980246	6.88e-05	0.0969047

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
33	33	mito	T	1	mito T	7938936	0.6465534	1
34	34	mito	T	2	mito T	12294342	0.8147126	1
47	47	mito	S	15	mito S	11142807	0.7268125	0
58	58	mito	T	26	mito T	6786437	0.4821588	1
59	59	mito	T	27	mito T	10802671	0.6646957	1
60	60	mito	T	28	mito T	9524968	0.5887804	1
64	64	mito	S	32	mito S	7138202	0.5016151	0

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

## final dispersion estimates

## fitting model and testing

##  chr "2 ENSRNOG00000059443.1 LOC103690028 protein_coding"
##  chr(0)

as.data.frame(de5[1:20,]) %>% kbl() %>% kable_paper("hover", full_width = F)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
2 ENSRNOG00000059443.1 LOC103690028 protein_coding	6.386214	30.0000000	4.1193787	7.282652	0.0000000	0.0000000
17 ENSRNOG00000056825.2 AABR07072539.4 processed_pseudogene	11.661670	6.5824326	1.5521786	4.240770	0.0000223	0.3578214
X ENSRNOG00000050700.2 AABR07037451.1 processed_pseudogene	1.519340	16.2028557	4.1292666	3.923906	0.0000871	0.7172777
1 ENSRNOG00000020025.4 LOC108348052 protein_coding	17.782509	3.2117296	0.8197475	3.917950	0.0000893	0.7172777
AABR07024203.1 ENSRNOG00000061777.1 AABR07024203.3 protein_coding	11.234939	-3.3373093	0.9087237	-3.672524	0.0002402	0.9999799
2 ENSRNOG00000020832.4 C2cd4d protein_coding	15.822890	2.8656038	0.8060464	3.555135	0.0003778	0.9999799
7 ENSRNOG00000050884.2 AABR07058745.1 protein_coding	50.195404	1.5965430	0.4942259	3.230391	0.0012362	0.9999799
3 ENSRNOG00000050864.2 Cpne1 protein_coding	14.072847	-2.5159879	0.7841478	-3.208563	0.0013340	0.9999799
18 ENSRNOG00000051569.1 LOC103690147 protein_coding	16.613644	2.5381822	0.7960416	3.188505	0.0014301	0.9999799
6 ENSRNOG00000051280.1 AABR07064502.1 processed_pseudogene	5.509108	5.4761178	1.7215016	3.181012	0.0014676	0.9999799
3 ENSRNOG00000007478.4 Cry2 protein_coding	471.596916	-0.5978197	0.1884367	-3.172523	0.0015112	0.9999799
4 ENSRNOG00000049639.2 LOC100909784 protein_coding	52.027683	1.1706273	0.3779391	3.097396	0.0019523	0.9999799
12 ENSRNOG00000057927.1 AABR07035796.1 lincRNA	58.135739	-0.9701573	0.3188153	-3.043008	0.0023423	0.9999799
KL568337.1 ENSRNOG00000057871.1 5_8S_rRNA rRNA	5.377655	-3.3118028	1.1068976	-2.991968	0.0027718	0.9999799
17 ENSRNOG00000052667.1 AABR07027555.1 processed_pseudogene	4.412753	5.1827233	1.7468898	2.966829	0.0030089	0.9999799
13 ENSRNOG00000056279.1 AC115369.1 lincRNA	199.904660	-0.5543759	0.1899197	-2.919001	0.0035116	0.9999799
13 ENSRNOG00000057562.1 U5 snRNA	9.045022	6.1759234	2.1248527	2.906518	0.0036548	0.9999799
14 ENSRNOG00000050406.1 AABR07014914.1 miRNA	102.689535	-0.7422696	0.2569182	-2.889128	0.0038631	0.9999799
13 ENSRNOG00000052639.1 LOC108352591 lincRNA	11.174558	2.3658472	0.8234583	2.873063	0.0040651	0.9999799
MT ENSRNOG00000029677.3 AY172581.6 Mt_tRNA	36.277712	1.5776326	0.5539265	2.848090	0.0043982	0.9999799

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	27	mean
MT ENSRNOG00000043866.1 AY172581.24 Mt_rRNA	65532.48480	103302.929	70376.072	143173.607	95596.273
MT ENSRNOG00000030478.3 AY172581.9 Mt_rRNA	68166.88500	94253.573	66097.117	127533.781	89012.839
10 ENSRNOG00000049695.2 Myh4 protein_coding	73162.30612	57759.766	36012.999	51313.467	54562.134
MT ENSRNOG00000034234.3 Mt-co1 protein_coding	27760.63254	40763.941	46883.092	57601.094	43252.190
19 ENSRNOG00000017786.6 Acta1 protein_coding	32867.41817	35902.178	28000.931	39461.381	34057.977
3 ENSRNOG00000006783.8 Neb protein_coding	15991.78610	16265.293	13783.228	16011.788	15513.024
15 ENSRNOG00000016983.5 Myh7 protein_coding	39.81012	9491.545	16050.108	22614.511	12048.993
MT ENSRNOG00000029707.3 Mt-nd4 protein_coding	9126.46772	11266.490	13322.592	13856.483	11893.008
1 ENSRNOG00000016837.4 Ckm protein_coding	12624.35786	12710.488	8756.212	13107.600	11799.664
8 ENSRNOG00000018184.8 Tpm1 protein_coding	14962.37709	9629.519	6935.752	5372.086	9224.933
10 ENSRNOG00000058068.1 Obscn protein_coding	8984.36928	10559.555	5807.829	10745.021	9024.193
1 ENSRNOG00000020332.8 Tnnt3 protein_coding	11001.61516	9555.417	7012.296	7849.811	8854.785
1 ENSRNOG00000047746.1 AABR07000398.1 protein_coding	7432.31378	10013.924	5021.958	12045.814	8628.502
MT ENSRNOG00000031979.3 Mt-atp6 protein_coding	5805.02841	8856.986	8647.247	9841.095	8287.589
MT ENSRNOG00000029971.3 Mt-nd5 protein_coding	6091.34969	8484.839	9423.286	7372.061	7842.884
MT ENSRNOG00000031033.4 Mt-nd2 protein_coding	6026.36481	7734.837	8826.612	7479.410	7516.806
1 ENSRNOG00000017645.6 Mylpf protein_coding	9366.22689	7831.735	6621.632	4611.283	7107.719
9 ENSRNOG00000013262.7 Myl1 protein_coding	8134.69807	6425.801	5028.650	4173.936	5940.771
MT ENSRNOG00000031766.3 Mt-cyb protein_coding	4330.41743	5386.463	6793.838	6764.630	5818.837
1 ENSRNOG00000019745.5 Actn3 protein_coding	9656.00605	6018.403	3966.229	2504.744	5536.345

# 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)

	20	26	28	32	7	mean
MT ENSRNOG00000043866.1 AY172581.24 Mt_rRNA	158053.630	158145.036	120971.027	88803.305	116980.109	128590.621
MT ENSRNOG00000030478.3 AY172581.9 Mt_rRNA	126811.529	117395.765	115558.939	137796.389	95135.763	118539.677
10 ENSRNOG00000049695.2 Myh4 protein_coding	67870.823	91464.506	66350.749	28310.810	64479.813	63695.340
MT ENSRNOG00000034234.3 Mt-co1 protein_coding	24143.686	20917.728	33563.379	19062.392	43784.838	28294.405
19 ENSRNOG00000017786.6 Acta1 protein_coding	25090.231	25867.914	45886.656	13388.493	20198.899	26086.439
3 ENSRNOG00000006783.8 Neb protein_coding	16435.901	22851.994	19473.485	24338.854	14385.281	19497.103
1 ENSRNOG00000047746.1 AABR07000398.1 protein_coding	12862.200	37968.805	30612.680	6824.503	6667.894	18987.217
1 ENSRNOG00000016837.4 Ckm protein_coding	8879.247	9639.926	15528.382	5646.266	14597.238	10858.212
15 ENSRNOG00000016983.5 Myh7 protein_coding	4828.164	9764.103	7522.180	15204.639	3675.789	8198.975
8 ENSRNOG00000018184.8 Tpm1 protein_coding	7759.914	8699.032	8612.078	5817.388	8226.563	7822.995
10 ENSRNOG00000058068.1 Obscn protein_coding	8718.745	8426.063	7970.619	5694.276	6489.291	7459.799
MT ENSRNOG00000029707.3 Mt-nd4 protein_coding	5900.972	6678.882	7573.295	7193.751	7924.565	7054.293
1 ENSRNOG00000020332.8 Tnnt3 protein_coding	6540.507	7796.425	7769.018	5800.355	7172.470	7015.755
MT ENSRNOG00000031979.3 Mt-atp6 protein_coding	5737.140	6467.459	7284.639	7344.178	7729.068	6912.497
1 ENSRNOG00000017645.6 Mylpf protein_coding	10198.437	4616.663	5993.887	4249.904	7677.033	6547.185
3 ENSRNOG00000015155.6 Tnnc2 protein_coding	5844.831	6097.463	5653.046	3236.298	8685.823	5903.492
MT ENSRNOG00000031766.3 Mt-cyb protein_coding	5046.097	4614.837	7223.435	3825.361	8751.170	5892.180
MT ENSRNOG00000029971.3 Mt-nd5 protein_coding	5929.410	4749.987	8432.616	4091.720	5821.454	5805.037
1 ENSRNOG00000052802.1 Aldoa protein_coding	8373.667	5556.689	6338.402	2967.876	5378.940	5723.115
1 ENSRNOG00000020557.7 Ryr1 protein_coding	6303.350	5543.083	5551.917	3461.411	4300.408	5032.034

# 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	11	12	13	14	16	17	18	19	21	22	23	24	25	29	3	30	31	4	5	6	8	9	mean
MT ENSRNOG00000043866.1 AY172581.24 Mt_rRNA	149740.771	140521.199	112912.776	90432.709	130645.108	109919.970	167770.260	225122.777	146316.179	179413.681	120425.159	171815.233	108694.330	155014.079	130115.959	92971.163	137715.899	127265.169	106936.933	158438.186	104784.387	88235.270	104604.318	133035.283
MT ENSRNOG00000030478.3 AY172581.9 Mt_rRNA	131108.710	149945.141	89966.596	109934.347	132868.836	113507.900	115621.269	180725.484	147278.915	217056.082	162290.953	167527.662	88908.063	139111.249	124143.158	105822.937	116400.986	149740.667	74681.961	201250.393	110581.446	78021.449	126470.807	131868.044
10 ENSRNOG00000049695.2 Myh4 protein_coding	41319.615	44973.305	67895.505	40608.684	50260.309	52639.233	53834.681	60329.612	56754.461	63505.595	53316.057	45731.297	71654.758	80910.627	37412.654	38575.221	78694.321	39132.047	75763.814	47704.054	38511.815	72470.527	36892.554	54299.598
MT ENSRNOG00000034234.3 Mt-co1 protein_coding	42494.672	45287.642	25971.541	34117.708	39616.493	33411.707	20474.443	34762.642	27268.106	42731.899	29603.385	35211.596	19038.883	32765.043	32074.744	24001.103	28595.124	33028.701	19483.939	69008.797	30447.769	16308.643	42109.344	32948.432
19 ENSRNOG00000017786.6 Acta1 protein_coding	32051.070	29717.123	28930.757	26112.701	23825.986	25834.074	22636.484	26729.431	22342.445	24228.131	25781.245	24531.816	24787.633	31491.556	19478.877	28308.117	23503.528	19805.499	28192.688	25923.501	16504.158	22091.308	21420.676	24966.470
3 ENSRNOG00000006783.8 Neb protein_coding	20751.164	18031.671	18943.388	19712.012	15765.207	26704.637	20689.644	28214.826	24040.888	20617.730	18882.818	17487.680	17539.411	27239.135	13873.188	10064.970	20209.093	13771.990	21286.314	18869.469	15835.261	16430.286	13570.219	19066.565
15 ENSRNOG00000016983.5 Myh7 protein_coding	24445.806	18335.812	7918.672	10920.299	13597.539	12539.089	16486.257	15396.703	14527.859	13243.229	3703.053	18619.346	5417.754	16000.234	13880.091	6869.012	8592.626	10965.078	3289.985	14969.633	14212.791	3009.103	8226.256	11963.749
1 ENSRNOG00000047746.1 AABR07000398.1 protein_coding	11618.883	7629.567	9939.861	12235.486	5870.149	7077.052	15272.867	9249.649	6408.034	23809.254	10323.266	13297.613	14777.998	18470.278	21826.798	6731.106	13482.958	10599.273	3456.447	9794.149	5725.101	7873.635	5117.162	10895.069
1 ENSRNOG00000016837.4 Ckm protein_coding	9041.923	12568.946	13148.389	13420.027	11932.222	9159.420	9510.592	8283.846	7978.927	12450.724	7716.590	6300.412	10685.825	10723.039	8913.875	11850.381	10704.491	7231.534	9951.961	10062.665	9145.936	11491.611	9498.986	10077.057
MT ENSRNOG00000029707.3 Mt-nd4 protein_coding	10990.765	9344.980	9702.389	8389.682	8226.199	8279.054	6386.713	8207.891	6976.830	10054.066	4869.040	11605.000	5104.697	9047.404	11144.789	5275.364	8030.187	7030.983	5970.828	12583.514	7888.611	3011.156	10713.377	8210.153
MT ENSRNOG00000031979.3 Mt-atp6 protein_coding	12122.001	8615.783	8472.782	10519.775	8923.397	7664.410	5616.194	8331.130	6586.636	10175.522	6049.870	6520.301	7021.557	8039.392	9234.867	7513.915	5598.808	9436.025	5221.714	11554.367	8660.772	4034.823	6852.102	7946.354
10 ENSRNOG00000058068.1 Obscn protein_coding	7410.582	4927.138	12244.612	8089.383	8556.778	9179.110	6204.813	9271.483	6427.235	8870.421	7082.295	7391.263	6173.368	10073.050	5204.206	4824.376	9546.306	4820.467	5416.171	5062.112	6353.665	7475.030	5786.130	7234.348
1 ENSRNOG00000020332.8 Tnnt3 protein_coding	6706.170	5567.826	8783.047	8904.134	6099.484	6189.682	5328.949	6804.493	6732.283	9325.355	7061.893	4688.227	10148.057	7067.174	5713.421	7318.939	5741.993	6299.777	7273.988	5829.653	6453.078	7044.795	4482.372	6763.687
MT ENSRNOG00000031766.3 Mt-cyb protein_coding	11344.638	6048.283	6079.296	8220.422	7666.163	4988.203	6358.528	6887.313	6445.615	9365.299	5338.689	6113.819	4182.381	7887.654	6124.670	5293.957	5620.574	4993.828	6536.386	9742.503	7336.217	4074.705	6100.529	6641.290
8 ENSRNOG00000018184.8 Tpm1 protein_coding	5337.411	4629.828	12277.784	6772.129	5853.576	6585.098	5665.120	7710.680	5355.808	8393.766	7977.525	4741.225	6264.043	7989.928	4174.731	8412.210	7282.033	4581.048	9887.128	3849.886	5580.562	7795.172	5001.830	6613.849
MT ENSRNOG00000029971.3 Mt-nd5 protein_coding	9213.737	7737.435	7415.710	7580.261	6820.254	10537.012	4945.853	7803.591	6731.391	6965.587	5806.536	5671.770	4454.312	7071.977	6189.793	3325.521	5352.576	5477.819	5227.341	10896.501	5861.881	3301.561	7063.666	6584.873
3 ENSRNOG00000015155.6 Tnnc2 protein_coding	5718.098	6577.086	9239.983	5230.299	5716.219	5117.589	4260.506	5446.792	4091.413	5804.313	3887.279	4520.201	5024.405	5157.058	4951.860	5699.687	6869.677	4675.299	8803.115	7515.961	3936.089	6279.263	4809.586	5623.121
1 ENSRNOG00000017645.6 Mylpf protein_coding	4738.342	4629.501	7863.004	6778.946	4861.412	5254.894	4195.354	4104.992	5680.912	6852.492	5124.485	4730.653	4942.195	6272.581	4439.091	7375.958	5660.905	5884.919	6740.231	5708.443	4914.775	6280.495	5796.469	5601.350
MT ENSRNOG00000030700.3 Mt-co3 protein_coding	5799.326	6067.711	5173.701	6703.132	7227.185	5193.084	3899.641	4896.239	4426.626	6389.612	4450.096	5533.308	2778.468	7720.524	5430.315	4710.764	4802.655	4832.243	4441.386	7383.763	4604.613	3070.121	5922.458	5280.738
1 ENSRNOG00000052802.1 Aldoa protein_coding	5002.804	3136.245	7381.938	5735.205	5457.096	5980.768	4314.311	4514.002	4377.757	6747.086	5150.095	4353.443	6459.005	5916.699	5348.211	4700.967	5941.405	4540.136	5262.108	5652.209	3883.178	4782.328	4819.744	5193.771

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)," ") , "[[", 2)
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 = 32852

## Note: no. genes in output = 18043

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

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] 71

head(res$enrichment_result,20) %>% kbl() %>% kable_paper("hover", full_width = F)

	set	setSize	pANOVA	s.dist	p.adjustANOVA
1440	rRNA processing in the nucleus and cytosol	166	0	-0.3243600	0e+00
668	Major pathway of rRNA processing in the nucleolus and cytosol	157	0	-0.3322112	0e+00
1110	SRP-dependent cotranslational protein targeting to membrane	89	0	-0.4181998	0e+00
860	Peptide chain elongation	66	0	-0.4797966	0e+00
1305	The citric acid (TCA) cycle and respiratory electron transport	165	0	0.3029130	0e+00
378	Eukaryotic Translation Termination	70	0	-0.4614159	0e+00
1074	Response of EIF2AK4 (GCN2) to amino acid deficiency	78	0	-0.4362978	0e+00
676	Metabolism	1943	0	0.0922760	0e+00
154	Cap-dependent Translation Initiation	93	0	-0.3984341	0e+00
377	Eukaryotic Translation Initiation	93	0	-0.3984341	0e+00
1130	Selenocysteine synthesis	69	0	-0.4616298	0e+00
1395	Viral mRNA Translation	66	0	-0.4674795	0e+00
376	Eukaryotic Translation Elongation	70	0	-0.4511815	0e+00
470	GTP hydrolysis and joining of the 60S ribosomal subunit	87	0	-0.3968257	0e+00
428	Formation of a pool of free 40S subunits	77	0	-0.4143765	0e+00
788	Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)	71	0	-0.4286543	0e+00
628	L13a-mediated translational silencing of Ceruloplasmin expression	85	0	-0.3914113	0e+00
787	Nonsense Mediated Decay (NMD) enhanced by the Exon Junction Complex (EJC)	90	0	-0.3579061	3e-07
789	Nonsense-Mediated Decay (NMD)	90	0	-0.3579061	3e-07
1438	rRNA processing	187	0	-0.2448791	6e-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] 44

nrow(subset(res$enrichment_result,p.adjustANOVA<0.05&s.dist>0))

## [1] 27

unlink("skeletal_mitch1.html")
mitch_report(res,outfile="skeletal_mitch1.html")

## Dataset saved as " /tmp/RtmperY9IA/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: /mnt/bfx6/bfx/adam_trewin/set2/mitch.knit.md

## /home/mdz/anaconda3/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/adam_trewin/set2/mitch.utf8.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmperY9IA/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/RtmperY9IA/rmarkdown-strc4115cec3ca5.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'

## 
## Output created: /tmp/RtmperY9IA/mitch_report.html

## [1] TRUE

Session info

sessionInfo()

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.3.4            pkgload_1.2.0              
##  [3] GGally_2.1.1                ggplot2_3.3.3              
##  [5] beeswarm_0.3.1              gtools_3.8.2               
##  [7] tibble_3.1.0                dplyr_1.0.5                
##  [9] echarts4r_0.4.0             gplots_3.1.1               
## [11] mitch_1.2.2                 DESeq2_1.30.0              
## [13] SummarizedExperiment_1.20.0 Biobase_2.50.0             
## [15] MatrixGenerics_1.2.0        matrixStats_0.58.0         
## [17] GenomicRanges_1.42.0        GenomeInfoDb_1.26.0        
## [19] IRanges_2.24.0              S4Vectors_0.28.0           
## [21] BiocGenerics_0.36.0         reshape2_1.4.4             
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-6           bit64_4.0.5            webshot_0.5.2         
##  [4] RColorBrewer_1.1-2     httr_1.4.2             rprojroot_2.0.2       
##  [7] tools_4.0.3            bslib_0.2.4            utf8_1.2.1            
## [10] R6_2.5.0               KernSmooth_2.23-18     DBI_1.1.1             
## [13] colorspace_2.0-0       withr_2.4.1            tidyselect_1.1.0      
## [16] gridExtra_2.3          bit_4.0.4              compiler_4.0.3        
## [19] rvest_1.0.0            xml2_1.3.2             desc_1.3.0            
## [22] DelayedArray_0.16.0    sass_0.3.1             caTools_1.18.1        
## [25] scales_1.1.1           genefilter_1.72.0      systemfonts_1.0.1     
## [28] stringr_1.4.0          digest_0.6.27          svglite_2.0.0         
## [31] rmarkdown_2.7          XVector_0.30.0         pkgconfig_2.0.3       
## [34] htmltools_0.5.1.1      highr_0.8              fastmap_1.1.0         
## [37] htmlwidgets_1.5.3      rlang_0.4.10           rstudioapi_0.13       
## [40] RSQLite_2.2.4          shiny_1.6.0            jquerylib_0.1.3       
## [43] generics_0.1.0         jsonlite_1.7.2         BiocParallel_1.24.1   
## [46] RCurl_1.98-1.3         magrittr_2.0.1         GenomeInfoDbData_1.2.4
## [49] Matrix_1.3-2           Rcpp_1.0.6             munsell_0.5.0         
## [52] fansi_0.4.2            lifecycle_1.0.0        stringi_1.5.3         
## [55] yaml_2.2.1             MASS_7.3-53.1          zlibbioc_1.36.0       
## [58] plyr_1.8.6             grid_4.0.3             blob_1.2.1            
## [61] promises_1.2.0.1       crayon_1.4.1           lattice_0.20-41       
## [64] splines_4.0.3          annotate_1.68.0        locfit_1.5-9.4        
## [67] knitr_1.31             pillar_1.5.1           geneplotter_1.68.0    
## [70] XML_3.99-0.6           glue_1.4.2             evaluate_0.14         
## [73] vctrs_0.3.6            httpuv_1.5.5           testthat_3.0.2        
## [76] gtable_0.3.0           purrr_0.3.4            reshape_0.8.8         
## [79] assertthat_0.2.1       cachem_1.0.4           xfun_0.22             
## [82] mime_0.10              xtable_1.8-4           later_1.1.0.1         
## [85] viridisLite_0.4.0      survival_3.2-10        AnnotationDbi_1.52.0  
## [88] memoise_2.0.0          ellipsis_0.3.1