Source: https://github.com/markziemann/mst1_reanalysis

Introduction

suppressPackageStartupMessages({
    library("reshape2")
    library("DESeq2")
    library("gplots")
    library("mitch")
    library("eulerr")
    library("kableExtra")
    library("beeswarm")
    library(ggplot2)
})

Import DE data

you <- read.table("Mst1TG_young_deseq.tsv",header=TRUE,sep="\t")

old <- read.table("Mst1TG_old_deseq.tsv",header=TRUE,sep="\t")

you_up <- rownames(subset(you,log2FoldChange > 0 & padj < 0.05))
you_dn <- rownames(subset(you,log2FoldChange < 0 & padj < 0.05))

old_up <- rownames(subset(old,log2FoldChange > 0 & padj < 0.05))
old_dn <- rownames(subset(old,log2FoldChange < 0 & padj < 0.05))

Euler diagram

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

v0 <- list("young up"=you_up, 
  "young dn"=you_dn,
  "old up"=old_up,
  "old dn"=old_dn)

plot(euler(v0),quantities = TRUE, edges = "gray", main="effect of Mst-1 overexpression")

Mitch

First fetch genesets.

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

Now set up the gene table gt.

gt <- as.data.frame(unique(rownames(old) ,rownames(you)))
gt$gene <- sapply(strsplit(gt[,1]," "),"[[",2)

Now execute mitch analysis.

m1 <- mitch_import(list("young"=you,"old"=old), DEtype="deseq2",geneTable=gt)
## Note: Mean no. genes in input = 17755
## Note: no. genes in output = 14877
## Note: estimated proportion of input genes in output = 0.838
mres1 <- mitch_calc(m1, genesets, priority="effect")
## Note: Enrichments with large effect sizes may not be 
##             statistically significant.
head(mres1$enrichment_result,20) %>% 
  kbl(caption = "Top gene pathway differences caused by Mst-1 overexpression") %>%
  kable_paper("hover", full_width = F)
Top gene pathway differences caused by Mst-1 overexpression
set setSize pMANOVA s.young s.old p.young p.old s.dist SD p.adjustMANOVA
165 COMPLEX I BIOGENESIS 55 0.0000000 -0.9369439 -0.9042259 0.0000000 0.0000000 1.3021092 0.0231351 0.0000000
322 FORMATION OF ATP BY CHEMIOSMOTIC COUPLING 18 0.0000000 -0.9281393 -0.8834003 0.0000000 0.0000000 1.2813425 0.0316353 0.0000000
149 CITRIC ACID CYCLE TCA CYCLE 22 0.0000000 -0.8959579 -0.8905541 0.0000000 0.0000000 1.2632606 0.0038210 0.0000000
816 RESPIRATORY ELECTRON TRANSPORT 101 0.0000000 -0.9086450 -0.8591642 0.0000000 0.0000000 1.2505194 0.0349882 0.0000000
183 CRISTAE FORMATION 31 0.0000000 -0.9061548 -0.8468361 0.0000000 0.0000000 1.2402612 0.0419446 0.0000000
817 RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CHEMIOSMOTIC COUPLING AND HEAT PRODUCTION BY UNCOUPLING PROTEINS 124 0.0000000 -0.8879965 -0.8342819 0.0000000 0.0000000 1.2184269 0.0379820 0.0000000
555 MITOCHONDRIAL FATTY ACID BETA OXIDATION 33 0.0000000 -0.8501466 -0.8532291 0.0000000 0.0000000 1.2044705 0.0021797 0.0000000
100 BRANCHED CHAIN AMINO ACID CATABOLISM 21 0.0000000 -0.8359105 -0.8450971 0.0000000 0.0000000 1.1886696 0.0064959 0.0000000
1043 THE CITRIC ACID TCA CYCLE AND RESPIRATORY ELECTRON TRANSPORT 171 0.0000000 -0.8413362 -0.8137348 0.0000000 0.0000000 1.1704746 0.0195171 0.0000000
558 MITOCHONDRIAL TRANSLATION 93 0.0000000 -0.8266769 -0.7564150 0.0000000 0.0000000 1.1205170 0.0496827 0.0000000
160 COLLAGEN CHAIN TRIMERIZATION 36 0.0000000 0.7396888 0.7579753 0.0000000 0.0000000 1.0590874 0.0129305 0.0000000
557 MITOCHONDRIAL PROTEIN IMPORT 63 0.0000000 -0.7830859 -0.7122820 0.0000000 0.0000000 1.0585694 0.0500659 0.0000000
748 PYRUVATE METABOLISM AND CITRIC ACID TCA CYCLE 51 0.0000000 -0.7273867 -0.7634257 0.0000000 0.0000000 1.0544715 0.0254834 0.0000000
385 GLYOXYLATE METABOLISM AND GLYCINE DEGRADATION 25 0.0000000 -0.7421492 -0.7394883 0.0000000 0.0000000 1.0476776 0.0018816 0.0000000
186 CROSSLINKING OF COLLAGEN FIBRILS 18 0.0000001 0.7541034 0.7000621 0.0000000 0.0000003 1.0289601 0.0382130 0.0000027
874 RRNA PROCESSING IN THE MITOCHONDRION 10 0.0001725 -0.6657967 -0.7546647 0.0002664 0.0000358 1.0063817 0.0628391 0.0018202
559 MITOCHONDRIAL TRNA AMINOACYLATION 21 0.0000000 -0.7061248 -0.6833731 0.0000000 0.0000001 0.9826551 0.0160879 0.0000012
493 KERATAN SULFATE DEGRADATION 11 0.0002225 0.6998276 0.6509424 0.0000583 0.0001851 0.9557639 0.0345671 0.0022025
269 ECM PROTEOGLYCANS 61 0.0000000 0.7080051 0.6314308 0.0000000 0.0000000 0.9486707 0.0541462 0.0000000
556 MITOCHONDRIAL IRON SULFUR CLUSTER BIOGENESIS 13 0.0000785 -0.6793181 -0.6408980 0.0000222 0.0000630 0.9339290 0.0271671 0.0009260

Scatterplot

res <- mres1$enrichment_result

sig <- subset(res,p.adjustMANOVA < 0.01)

plot(res$s.young,res$s.old,bty="none",
  xlim=c(-1,1), ylim=c(-1,1),
  xlab="young",ylab="old",
  main="enrichment score, effect of Mst-1 overexpression") 

abline(v=0,lty=2)
abline(h=0,lty=2)

points(sig$s.young, sig$s.old, pch=19,col="black")

message("downregulated gene sets")
## downregulated gene sets
dn <- subset(sig,s.young < -0.5 & s.old < -0.5)$set
dn <- unname(sapply(X=dn,function(x) {substring(x, first=1, last=50) } ))
dn
##  [1] "COMPLEX I BIOGENESIS"                              
##  [2] "FORMATION OF ATP BY CHEMIOSMOTIC COUPLING"         
##  [3] "CITRIC ACID CYCLE TCA CYCLE"                       
##  [4] "RESPIRATORY ELECTRON TRANSPORT"                    
##  [5] "CRISTAE FORMATION"                                 
##  [6] "RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CH"
##  [7] "MITOCHONDRIAL FATTY ACID BETA OXIDATION"           
##  [8] "BRANCHED CHAIN AMINO ACID CATABOLISM"              
##  [9] "THE CITRIC ACID TCA CYCLE AND RESPIRATORY ELECTRON"
## [10] "MITOCHONDRIAL TRANSLATION"                         
## [11] "MITOCHONDRIAL PROTEIN IMPORT"                      
## [12] "PYRUVATE METABOLISM AND CITRIC ACID TCA CYCLE"     
## [13] "GLYOXYLATE METABOLISM AND GLYCINE DEGRADATION"     
## [14] "RRNA PROCESSING IN THE MITOCHONDRION"              
## [15] "MITOCHONDRIAL TRNA AMINOACYLATION"                 
## [16] "MITOCHONDRIAL IRON SULFUR CLUSTER BIOGENESIS"      
## [17] "REGULATION OF PYRUVATE DEHYDROGENASE PDH COMPLEX"  
## [18] "METABOLISM OF COFACTORS"                           
## [19] "TRIGLYCERIDE CATABOLISM"                           
## [20] "PYRUVATE METABOLISM"                               
## [21] "BETA OXIDATION OF VERY LONG CHAIN FATTY ACIDS"     
## [22] "DEGRADATION OF CYSTEINE AND HOMOCYSTEINE"          
## [23] "DEFECTS IN COBALAMIN B12 METABOLISM"               
## [24] "MITOCHONDRIAL BIOGENESIS"                          
## [25] "PROTEIN LOCALIZATION"                              
## [26] "HEME BIOSYNTHESIS"                                 
## [27] "MITOCHONDRIAL CALCIUM ION TRANSPORT"
message("upregulated gene sets")
## upregulated gene sets
up <- subset(sig,s.young > 0.5 & s.old > 0.5)$set
up <- unname(sapply(X=up,function(x) {substring(x, first=1, last=45) } ))
up
##  [1] "COLLAGEN CHAIN TRIMERIZATION"                 
##  [2] "CROSSLINKING OF COLLAGEN FIBRILS"             
##  [3] "KERATAN SULFATE DEGRADATION"                  
##  [4] "ECM PROTEOGLYCANS"                            
##  [5] "COLLAGEN BIOSYNTHESIS AND MODIFYING ENZYMES"  
##  [6] "CS DS DEGRADATION"                            
##  [7] "MET ACTIVATES PTK2 SIGNALING"                 
##  [8] "CHONDROITIN SULFATE BIOSYNTHESIS"             
##  [9] "ANCHORING FIBRIL FORMATION"                   
## [10] "SYNDECAN INTERACTIONS"                        
## [11] "LAMININ INTERACTIONS"                         
## [12] "ELASTIC FIBRE FORMATION"                      
## [13] "ASSEMBLY OF COLLAGEN FIBRILS AND OTHER MULTIM"
## [14] "COLLAGEN DEGRADATION"                         
## [15] "COLLAGEN FORMATION"                           
## [16] "NON INTEGRIN MEMBRANE ECM INTERACTIONS"       
## [17] "MOLECULES ASSOCIATED WITH ELASTIC FIBRES"
points(-0.75,-0.75,cex=30,col="red")

pos <- data.frame(-seq(0.35,1,0.05))
pos$x <- 0.05

text(0.35,-0.3,"DOWNREGULATED",cex=1,col="red")
text(pos[,2],pos[,1], head(dn,  nrow(pos)) , pos=4,cex=0.8)

points(0.75,0.75,cex=30,col="blue")

pos <- data.frame(seq(0.9,0.25,-0.05))
pos$x <- -1

text(-0.75,0.95,"UPREGULATED",cex=1,col="blue")
text(pos[,2],pos[,1], head(up,  nrow(pos)) , pos=4,cex=0.8)

Get counts for heatmaps

Need to read these in from the kallisto output, clean up the rownames and then cleanup the column names.

Then combine the 2 tables into a single one, then calculate RPM.

tmp <- read.table("3col_young.tsv.gz",header=F)
x <- as.matrix(acast(tmp, V2~V1, value.var="V3", fun.aggregate = sum))
x <- as.data.frame(x)
accession <- sapply((strsplit(rownames(x),"\\|")),"[[",2)
symbol<-sapply((strsplit(rownames(x),"\\|")),"[[",6)
x$geneid <- paste(accession,symbol)
xy <- aggregate(. ~ geneid,x,sum)
rownames(xy) <- xy$geneid
xy$geneid = NULL
colnames(xy) <- gsub("D","Mst1tg",colnames(xy))
colnames(xy) <- gsub("E","wt",colnames(xy))
colnames(xy) <- paste(colnames(xy) ,"_young",sep="")
rpmy <- xy / colSums(xy) * 1000000
rpmy2 <- rpmy
rpmy <- rpmy[which(rowMeans(xy) > 10),]

tmp <- read.table("3col_old.tsv.gz",header=F)
x <- as.matrix(acast(tmp, V2~V1, value.var="V3", fun.aggregate = sum))
x <- as.data.frame(x)
accession <- sapply((strsplit(rownames(x),"\\|")),"[[",2)
symbol<-sapply((strsplit(rownames(x),"\\|")),"[[",6)
x$geneid <- paste(accession,symbol)
xo <- aggregate(. ~ geneid,x,sum)
rownames(xo) <- xo$geneid
xo$geneid = NULL
xo <-  xo[,grep("Gal",colnames(xo),invert=TRUE)] 
colnames(xo) <- sapply(strsplit(colnames(xo),"\\."),"[[",1)
colnames(xo) <- paste(colnames(xo) ,"_old",sep="")
rpmo <- xo / colSums(xo) * 1000000
rpmo2 <- rpmo
rpmo <- rpmo[which(rowMeans(xo) > 10),]

xx <- cbind(xy,xo)
rpm <- xx / colSums(xx) * 1000000 
rpm2 <- rpm
head(rpm)
##                             Mst1tg1_young Mst1tg2_young Mst1tg3_young
## ENSMUSG00000000001.4 Gnai3     65.0401759   48.20020933    52.3868676
## ENSMUSG00000000003.15 Pbsn      0.0000000    0.00000000     0.0000000
## ENSMUSG00000000028.15 Cdc45     4.2879384    2.24580409     2.2252471
## ENSMUSG00000000031.16 H19     225.2088355  109.01299021   228.6548749
## ENSMUSG00000000037.17 Scml2     2.0160082    1.07277773     0.8044967
## ENSMUSG00000000049.11 Apoh      0.7362421    0.03860601     0.2236595
##                             Mst1tg4_young  wt1_young  wt2_young   wt3_young
## ENSMUSG00000000001.4 Gnai3     59.9427146 30.1402083  29.521231  43.9295654
## ENSMUSG00000000003.15 Pbsn      0.0000000  0.0000000   0.000000   0.0000000
## ENSMUSG00000000028.15 Cdc45     3.0369973  2.3935720   2.012936   2.5077732
## ENSMUSG00000000031.16 H19     210.9220024 95.5398937 134.777873 117.9208203
## ENSMUSG00000000037.17 Scml2     1.5098869  0.6036634   1.290992   0.6114870
## ENSMUSG00000000049.11 Apoh      0.1700184  0.7291144   1.190612   0.1930003
##                               wt4_young MstTg1_old MstTg2_old MstTg3_old
## ENSMUSG00000000001.4 Gnai3   31.2522740 27.8046569  24.264894 24.8642023
## ENSMUSG00000000003.15 Pbsn    0.0000000  0.0000000   0.000000  0.0000000
## ENSMUSG00000000028.15 Cdc45   1.6769627  3.2529203   1.817730  2.3549100
## ENSMUSG00000000031.16 H19   121.2129335  5.7540424   4.968648  7.2702675
## ENSMUSG00000000037.17 Scml2   0.4944992  1.4264694   1.467030  1.7445431
## ENSMUSG00000000049.11 Apoh    0.5034071  0.7307283   0.259355  0.2799127
##                             MstTg4_old MstTg5_old MstTg6_old MstTg7_old
## ENSMUSG00000000001.4 Gnai3   30.240133 21.1450025 16.5070827 33.5477577
## ENSMUSG00000000003.15 Pbsn    0.000000  0.0000000  0.0000000  0.0000000
## ENSMUSG00000000028.15 Cdc45   2.758618  2.5623960  2.5310887  3.5470401
## ENSMUSG00000000031.16 H19     7.785368  7.6128741  8.8525063  5.0357753
## ENSMUSG00000000037.17 Scml2   1.785192  1.5311384  1.6027473  1.8978342
## ENSMUSG00000000049.11 Apoh    0.520534  0.4220775  0.1260006  0.3105408
##                                wt1_old    wt2_old    wt3_old    wt4_old
## ENSMUSG00000000001.4 Gnai3  14.6556378 19.2265462 15.3591129 15.5355213
## ENSMUSG00000000003.15 Pbsn   0.0000000  0.0000000  0.0000000  0.0000000
## ENSMUSG00000000028.15 Cdc45  1.7045155  1.5164998  0.8228980  1.1110777
## ENSMUSG00000000031.16 H19    4.7801601 16.2893775  5.5206605  5.5020360
## ENSMUSG00000000037.17 Scml2  0.5393648  0.7307279  0.6483869  0.8397385
## ENSMUSG00000000049.11 Apoh   0.2085732  0.2831035  0.2012216  0.1383205
##                                wt5_old    wt6_old
## ENSMUSG00000000001.4 Gnai3  19.2343253 22.8547590
## ENSMUSG00000000003.15 Pbsn   0.0000000  0.0000000
## ENSMUSG00000000028.15 Cdc45  1.8476114  1.6264769
## ENSMUSG00000000031.16 H19    7.3532953  3.3539253
## ENSMUSG00000000037.17 Scml2  0.6406568  0.6523002
## ENSMUSG00000000049.11 Apoh   0.0000000  0.2884580
dim(rpm)
## [1] 55385    21
rpm <- rpm[which(rowMeans(xx) > 10),]
dim(rpm)
## [1] 19083    21

Mito encoded genes

mty <- rpmy2[grep(" mt-",rownames(rpmy2)),]
mty <- mty[grep("mt-T",rownames(mty),invert=TRUE),]
mtywt <- mty[,grep("wt",colnames(mty))]
mtytg <- mty[,grep("tg",colnames(mty))]

mto <- rpmo2[grep(" mt-",rownames(rpmo2)),]
mto <- mto[grep("mt-T",rownames(mto),invert=TRUE),]
mtowt <- mto[,grep("wt",colnames(mto))]
mtotg <- mto[,grep("Tg",colnames(mto))]

myboxplot <- function(GNAME) {

dat <- list("young wt"=as.numeric(mtywt[GNAME,]),
    "young tg"=as.numeric(mtytg[GNAME,]) ,
    "old wt"=as.numeric(mtowt[GNAME,]),
    "old tg"=as.numeric(mtotg[GNAME,]) )

boxplot(dat, main=GNAME, col="white", 
  log="y", ylab="RPM (log scale)")

beeswarm(dat,add=TRUE,pch=19,cex=1.2)
}

GNAMES <- intersect( rownames(mto) , rownames(mty) )

myboxplot(GNAMES[1])

par(mfrow=c(2,2))

sapply(GNAMES,myboxplot)

##        ENSMUSG00000064337.1 mt-Rnr1 ENSMUSG00000064339.1 mt-Rnr2
## x      numeric,21                   numeric,21                  
## y      numeric,21                   numeric,21                  
## pch    numeric,21                   numeric,21                  
## col    character,21                 character,21                
## bg     logical,21                   logical,21                  
## cex    numeric,21                   numeric,21                  
## x.orig character,21                 character,21                
## y.orig numeric,21                   numeric,21                  
##        ENSMUSG00000064341.1 mt-Nd1 ENSMUSG00000064345.1 mt-Nd2
## x      numeric,21                  numeric,21                 
## y      numeric,21                  numeric,21                 
## pch    numeric,21                  numeric,21                 
## col    character,21                character,21               
## bg     logical,21                  logical,21                 
## cex    numeric,21                  numeric,21                 
## x.orig character,21                character,21               
## y.orig numeric,21                  numeric,21                 
##        ENSMUSG00000064351.1 mt-Co1 ENSMUSG00000064354.1 mt-Co2
## x      numeric,21                  numeric,21                 
## y      numeric,21                  numeric,21                 
## pch    numeric,21                  numeric,21                 
## col    character,21                character,21               
## bg     logical,21                  logical,21                 
## cex    numeric,21                  numeric,21                 
## x.orig character,21                character,21               
## y.orig numeric,21                  numeric,21                 
##        ENSMUSG00000064356.3 mt-Atp8 ENSMUSG00000064357.1 mt-Atp6
## x      numeric,21                   numeric,21                  
## y      numeric,21                   numeric,21                  
## pch    numeric,21                   numeric,21                  
## col    character,21                 character,21                
## bg     logical,21                   logical,21                  
## cex    numeric,21                   numeric,21                  
## x.orig character,21                 character,21                
## y.orig numeric,21                   numeric,21                  
##        ENSMUSG00000064358.1 mt-Co3 ENSMUSG00000064360.1 mt-Nd3
## x      numeric,21                  numeric,21                 
## y      numeric,21                  numeric,21                 
## pch    numeric,21                  numeric,21                 
## col    character,21                character,21               
## bg     logical,21                  logical,21                 
## cex    numeric,21                  numeric,21                 
## x.orig character,21                character,21               
## y.orig numeric,21                  numeric,21                 
##        ENSMUSG00000064363.1 mt-Nd4 ENSMUSG00000064367.1 mt-Nd5
## x      numeric,21                  numeric,21                 
## y      numeric,21                  numeric,21                 
## pch    numeric,21                  numeric,21                 
## col    character,21                character,21               
## bg     logical,21                  logical,21                 
## cex    numeric,21                  numeric,21                 
## x.orig character,21                character,21               
## y.orig numeric,21                  numeric,21                 
##        ENSMUSG00000064368.1 mt-Nd6 ENSMUSG00000064370.1 mt-Cytb
## x      numeric,21                  numeric,21                  
## y      numeric,21                  numeric,21                  
## pch    numeric,21                  numeric,21                  
## col    character,21                character,21                
## bg     logical,21                  logical,21                  
## cex    numeric,21                  numeric,21                  
## x.orig character,21                character,21                
## y.orig numeric,21                  numeric,21                  
##        ENSMUSG00000065947.3 mt-Nd4l
## x      numeric,21                  
## y      numeric,21                  
## pch    numeric,21                  
## col    character,21                
## bg     logical,21                  
## cex    numeric,21                  
## x.orig character,21                
## y.orig numeric,21

Try with ggplot.

mty <- rpmy2[grep(" mt-",rownames(rpmy2)),]
mtywt <- mty[,grep("wt",colnames(mty))]
mtywtl <- melt(mtywt)
## No id variables; using all as measure variables
mtywtl$gene <- rep(rownames(mtywt),ncol(mtywt))
mtywtl <- mtywtl[,c(3,1,2)]
colnames(mtywtl) <- c("gene","sample","quant")
mtywtl$sample <- "wt"

mtytg <- mty[,grep("tg",colnames(mty))]
mtytgl <- melt(mtytg)
## No id variables; using all as measure variables
mtytgl$gene <- rep(rownames(mtytg),ncol(mtytg))
mtytgl <- mtytgl[,c(3,1,2)]
colnames(mtytgl) <- c("gene","sample","quant")
mtytgl$sample <- "tg"

mtl <- rbind(mtywtl, mtytgl)

ggplot(mtl, aes(x = factor(gene), y = quant, fill = factor(sample))) + 
  geom_boxplot() 

mtl <- mtl[grep("mt-T",mtl$gene,invert=TRUE),]

ggplot(mtl, aes(x = factor(gene), y = quant, fill = factor(sample))) +
  geom_boxplot()

Need to show the RQ only.

mty <- rpmy2[grep(" mt-",rownames(rpmy2)),]
mty <- mty[grep("mt-T",rownames(mty),invert=TRUE),]
mtywt <- mty[,grep("wt",colnames(mty))]
mtywt_means <- rowMeans(mtywt)
mtytg <- mty[,grep("tg",colnames(mty))]
mtytg_means <- rowMeans(mtytg)
mtywt <- mtywt / mtywt_means
mtytg <- mtytg / mtywt_means

mtywtl <- melt(mtywt)
## No id variables; using all as measure variables
mtywtl$gene <- rep(rownames(mtywt),ncol(mtywt))
mtywtl <- mtywtl[,c(3,1,2)]
colnames(mtywtl) <- c("gene","sample","quant")
mtywtl$sample <- "wt"

mtytgl <- melt(mtytg)
## No id variables; using all as measure variables
mtytgl$gene <- rep(rownames(mtytg),ncol(mtytg))
mtytgl <- mtytgl[,c(3,1,2)]
colnames(mtytgl) <- c("gene","sample","quant")
mtytgl$sample <- "tg"

mtl <- rbind(mtywtl, mtytgl)

mtl$gene <- sapply(strsplit(mtl$gene , " "),"[[",2)

ggplot(mtl, aes(x = factor(gene), y = quant, fill = factor(sample))) +
  coord_flip() +
  geom_boxplot() +
  ggtitle("young") +
  theme_bw() +
  theme(text = element_text(size = 20))  

mto <- rpmo2[grep(" mt-",rownames(rpmo2)),]
mto <- mto[grep("mt-T",rownames(mto),invert=TRUE),]
mtowt <- mto[,grep("wt",colnames(mto))]
mtowt_means <- rowMeans(mtowt)
mtotg <- mto[,grep("Tg",colnames(mto))]
mtotg_means <- rowMeans(mtotg)
mtowt <- mtowt / mtowt_means
mtotg <- mtotg / mtowt_means

mtowtl <- melt(mtowt)
## No id variables; using all as measure variables
mtowtl$gene <- rep(rownames(mtowt),ncol(mtowt))
mtowtl <- mtowtl[,c(3,1,2)]
colnames(mtowtl) <- c("gene","sample","quant")
mtowtl$sample <- "wt"

mtotgl <- melt(mtotg)
## No id variables; using all as measure variables
mtotgl$gene <- rep(rownames(mtotg),ncol(mtotg))
mtotgl <- mtotgl[,c(3,1,2)]
colnames(mtotgl) <- c("gene","sample","quant")
mtotgl$sample <- "tg"

mtl <- rbind(mtowtl, mtotgl)

mtl$gene <- sapply(strsplit(mtl$gene , " "),"[[",2)

ggplot(mtl, aes(x = factor(gene), y = quant, fill = factor(sample))) +
  coord_flip() +
  geom_boxplot() +
  ggtitle("old") +
  theme_bw() +
  theme(text = element_text(size = 20))

Heatmaps

Here, I’m going to create heatmaps for some particular gene sets.

  • Respiratory electron transport

  • Collagen biosynthesis

  • mtDNA repair

  • Mitochondrial RNA binding

They may look a bit messy because combining separate experiments can be problematic.

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

rpm_genes <- sapply(strsplit(rownames(rpm)," "),"[[",2)

SETS <- unique( c( names(genesets)[grep("RESPIRATORY",names(genesets))],
    names(genesets)[grep("COLLAGEN",names(genesets))], 
    names(genesets)[grep("REPAIR",names(genesets))],
    names(genesets)[grep("MITOCHONDRIAL",names(genesets))] ) )

myheatmap <- function(SET) {
    message(SET)
    geneset_members <- genesets[[which(names(genesets)==SET)]]

    if ( length(which(rpm_genes %in% geneset_members)) > 10 ) {

        rpm_subset <- rpm[which(rpm_genes %in% geneset_members),]

        gnames <- sapply(strsplit(rownames(rpm_subset)," "),"[[",2)
        gnames[which(duplicated(gnames) )] <- paste( gnames[which(duplicated(gnames) )] , "2" )
        rownames(rpm_subset) <- gnames
        heatmap.2(as.matrix(rpm_subset), trace="none", scale="row", 
            col=colfunc(25) , dendrogram="none", margins = c(10,15))
        mtext(SET)
    }
}

sapply(SETS,myheatmap)
## RESPIRATORY ELECTRON TRANSPORT
## RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CHEMIOSMOTIC COUPLING AND HEAT PRODUCTION BY UNCOUPLING PROTEINS

## THE CITRIC ACID TCA CYCLE AND RESPIRATORY ELECTRON TRANSPORT

## ASSEMBLY OF COLLAGEN FIBRILS AND OTHER MULTIMERIC STRUCTURES

## COLLAGEN BIOSYNTHESIS AND MODIFYING ENZYMES

## COLLAGEN CHAIN TRIMERIZATION

## COLLAGEN DEGRADATION

## COLLAGEN FORMATION

## CROSSLINKING OF COLLAGEN FIBRILS

## PLATELET ADHESION TO EXPOSED COLLAGEN

## BASE EXCISION REPAIR

## BASE EXCISION REPAIR AP SITE FORMATION

## DISEASES OF BASE EXCISION REPAIR
## DISEASES OF DNA REPAIR

## DISEASES OF MISMATCH REPAIR MMR
## DNA DOUBLE STRAND BREAK REPAIR

## DNA REPAIR

## GAP FILLING DNA REPAIR SYNTHESIS AND LIGATION IN GG NER

## GLOBAL GENOME NUCLEOTIDE EXCISION REPAIR GG NER

## HOMOLOGY DIRECTED REPAIR

## MISMATCH REPAIR

## NUCLEOTIDE EXCISION REPAIR

## PCNA DEPENDENT LONG PATCH BASE EXCISION REPAIR

## POLB DEPENDENT LONG PATCH BASE EXCISION REPAIR
## PROTEIN REPAIR
## SUMOYLATION OF DNA DAMAGE RESPONSE AND REPAIR PROTEINS

## TP53 REGULATES TRANSCRIPTION OF DNA REPAIR GENES

## TRANSCRIPTION COUPLED NUCLEOTIDE EXCISION REPAIR TC NER

## MITOCHONDRIAL BIOGENESIS

## MITOCHONDRIAL CALCIUM ION TRANSPORT

## MITOCHONDRIAL FATTY ACID BETA OXIDATION

## MITOCHONDRIAL FATTY ACID BETA OXIDATION OF SATURATED FATTY ACIDS
## MITOCHONDRIAL FATTY ACID BETA OXIDATION OF UNSATURATED FATTY ACIDS
## MITOCHONDRIAL IRON SULFUR CLUSTER BIOGENESIS

## MITOCHONDRIAL PROTEIN IMPORT

## MITOCHONDRIAL TRANSLATION

## MITOCHONDRIAL TRNA AMINOACYLATION

## MITOCHONDRIAL UNCOUPLING
## TRANSCRIPTIONAL ACTIVATION OF MITOCHONDRIAL BIOGENESIS

## $`RESPIRATORY ELECTRON TRANSPORT`
## NULL
## 
## $`RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CHEMIOSMOTIC COUPLING AND HEAT PRODUCTION BY UNCOUPLING PROTEINS`
## NULL
## 
## $`THE CITRIC ACID TCA CYCLE AND RESPIRATORY ELECTRON TRANSPORT`
## NULL
## 
## $`ASSEMBLY OF COLLAGEN FIBRILS AND OTHER MULTIMERIC STRUCTURES`
## NULL
## 
## $`COLLAGEN BIOSYNTHESIS AND MODIFYING ENZYMES`
## NULL
## 
## $`COLLAGEN CHAIN TRIMERIZATION`
## NULL
## 
## $`COLLAGEN DEGRADATION`
## NULL
## 
## $`COLLAGEN FORMATION`
## NULL
## 
## $`CROSSLINKING OF COLLAGEN FIBRILS`
## NULL
## 
## $`PLATELET ADHESION TO EXPOSED COLLAGEN`
## NULL
## 
## $`BASE EXCISION REPAIR`
## NULL
## 
## $`BASE EXCISION REPAIR AP SITE FORMATION`
## NULL
## 
## $`DISEASES OF BASE EXCISION REPAIR`
## NULL
## 
## $`DISEASES OF DNA REPAIR`
## NULL
## 
## $`DISEASES OF MISMATCH REPAIR MMR`
## NULL
## 
## $`DNA DOUBLE STRAND BREAK REPAIR`
## NULL
## 
## $`DNA REPAIR`
## NULL
## 
## $`GAP FILLING DNA REPAIR SYNTHESIS AND LIGATION IN GG NER`
## NULL
## 
## $`GLOBAL GENOME NUCLEOTIDE EXCISION REPAIR GG NER`
## NULL
## 
## $`HOMOLOGY DIRECTED REPAIR`
## NULL
## 
## $`MISMATCH REPAIR`
## NULL
## 
## $`NUCLEOTIDE EXCISION REPAIR`
## NULL
## 
## $`PCNA DEPENDENT LONG PATCH BASE EXCISION REPAIR`
## NULL
## 
## $`POLB DEPENDENT LONG PATCH BASE EXCISION REPAIR`
## NULL
## 
## $`PROTEIN REPAIR`
## NULL
## 
## $`SUMOYLATION OF DNA DAMAGE RESPONSE AND REPAIR PROTEINS`
## NULL
## 
## $`TP53 REGULATES TRANSCRIPTION OF DNA REPAIR GENES`
## NULL
## 
## $`TRANSCRIPTION COUPLED NUCLEOTIDE EXCISION REPAIR TC NER`
## NULL
## 
## $`MITOCHONDRIAL BIOGENESIS`
## NULL
## 
## $`MITOCHONDRIAL CALCIUM ION TRANSPORT`
## NULL
## 
## $`MITOCHONDRIAL FATTY ACID BETA OXIDATION`
## NULL
## 
## $`MITOCHONDRIAL FATTY ACID BETA OXIDATION OF SATURATED FATTY ACIDS`
## NULL
## 
## $`MITOCHONDRIAL FATTY ACID BETA OXIDATION OF UNSATURATED FATTY ACIDS`
## NULL
## 
## $`MITOCHONDRIAL IRON SULFUR CLUSTER BIOGENESIS`
## NULL
## 
## $`MITOCHONDRIAL PROTEIN IMPORT`
## NULL
## 
## $`MITOCHONDRIAL TRANSLATION`
## NULL
## 
## $`MITOCHONDRIAL TRNA AMINOACYLATION`
## NULL
## 
## $`MITOCHONDRIAL UNCOUPLING`
## NULL
## 
## $`TRANSCRIPTIONAL ACTIVATION OF MITOCHONDRIAL BIOGENESIS`
## NULL

They may look a bit messy because combining separate experiments can be problematic.

rpmy_genes <- sapply(strsplit(rownames(rpmy)," "),"[[",2)
rpmo_genes <- sapply(strsplit(rownames(rpmo)," "),"[[",2)

myheatmap2 <- function(SET) {
    message(SET)
    geneset_members <- genesets[[which(names(genesets)==SET)]]
    if ( length(which(rpmy_genes %in% geneset_members)) > 10 ) {
        rpm_subset <- rpmy[which(rpmy_genes %in% geneset_members),]
        gnames <- sapply(strsplit(rownames(rpm_subset)," "),"[[",2)
        gnames[which(duplicated(gnames) )] <- paste( gnames[which(duplicated(gnames) )] , "2" )
        rownames(rpm_subset) <- gnames
        heatmap.2(as.matrix(rpm_subset), trace="none", scale="row",
            col=colfunc(25) , dendrogram="none", margins = c(10,18) )
        mtext(SET)
    }
    if ( length(which(rpmo_genes %in% geneset_members)) > 10 ) {
        rpm_subset <- rpmo[which(rpmo_genes %in% geneset_members),]
        gnames <- sapply(strsplit(rownames(rpm_subset)," "),"[[",2)
        gnames[which(duplicated(gnames) )] <- paste( gnames[which(duplicated(gnames) )] , "2" )
        rownames(rpm_subset) <- gnames
        heatmap.2(as.matrix(rpm_subset), trace="none", scale="row",
            col=colfunc(25) , dendrogram="none", margins = c(10,18))
        mtext(SET)
    }
}

sapply(SETS,myheatmap2)
## RESPIRATORY ELECTRON TRANSPORT

## RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CHEMIOSMOTIC COUPLING AND HEAT PRODUCTION BY UNCOUPLING PROTEINS

## THE CITRIC ACID TCA CYCLE AND RESPIRATORY ELECTRON TRANSPORT

## ASSEMBLY OF COLLAGEN FIBRILS AND OTHER MULTIMERIC STRUCTURES

## COLLAGEN BIOSYNTHESIS AND MODIFYING ENZYMES

## COLLAGEN CHAIN TRIMERIZATION

## COLLAGEN DEGRADATION

## COLLAGEN FORMATION

## CROSSLINKING OF COLLAGEN FIBRILS

## PLATELET ADHESION TO EXPOSED COLLAGEN

## BASE EXCISION REPAIR

## BASE EXCISION REPAIR AP SITE FORMATION

## DISEASES OF BASE EXCISION REPAIR
## DISEASES OF DNA REPAIR

## DISEASES OF MISMATCH REPAIR MMR
## DNA DOUBLE STRAND BREAK REPAIR

## DNA REPAIR

## GAP FILLING DNA REPAIR SYNTHESIS AND LIGATION IN GG NER

## GLOBAL GENOME NUCLEOTIDE EXCISION REPAIR GG NER

## HOMOLOGY DIRECTED REPAIR

## MISMATCH REPAIR

## NUCLEOTIDE EXCISION REPAIR

## PCNA DEPENDENT LONG PATCH BASE EXCISION REPAIR

## POLB DEPENDENT LONG PATCH BASE EXCISION REPAIR
## PROTEIN REPAIR
## SUMOYLATION OF DNA DAMAGE RESPONSE AND REPAIR PROTEINS

## TP53 REGULATES TRANSCRIPTION OF DNA REPAIR GENES

## TRANSCRIPTION COUPLED NUCLEOTIDE EXCISION REPAIR TC NER

## MITOCHONDRIAL BIOGENESIS

## MITOCHONDRIAL CALCIUM ION TRANSPORT

## MITOCHONDRIAL FATTY ACID BETA OXIDATION

## MITOCHONDRIAL FATTY ACID BETA OXIDATION OF SATURATED FATTY ACIDS
## MITOCHONDRIAL FATTY ACID BETA OXIDATION OF UNSATURATED FATTY ACIDS
## MITOCHONDRIAL IRON SULFUR CLUSTER BIOGENESIS

## MITOCHONDRIAL PROTEIN IMPORT

## MITOCHONDRIAL TRANSLATION

## MITOCHONDRIAL TRNA AMINOACYLATION

## MITOCHONDRIAL UNCOUPLING
## TRANSCRIPTIONAL ACTIVATION OF MITOCHONDRIAL BIOGENESIS

## $`RESPIRATORY ELECTRON TRANSPORT`
## NULL
## 
## $`RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS BY CHEMIOSMOTIC COUPLING AND HEAT PRODUCTION BY UNCOUPLING PROTEINS`
## NULL
## 
## $`THE CITRIC ACID TCA CYCLE AND RESPIRATORY ELECTRON TRANSPORT`
## NULL
## 
## $`ASSEMBLY OF COLLAGEN FIBRILS AND OTHER MULTIMERIC STRUCTURES`
## NULL
## 
## $`COLLAGEN BIOSYNTHESIS AND MODIFYING ENZYMES`
## NULL
## 
## $`COLLAGEN CHAIN TRIMERIZATION`
## NULL
## 
## $`COLLAGEN DEGRADATION`
## NULL
## 
## $`COLLAGEN FORMATION`
## NULL
## 
## $`CROSSLINKING OF COLLAGEN FIBRILS`
## NULL
## 
## $`PLATELET ADHESION TO EXPOSED COLLAGEN`
## NULL
## 
## $`BASE EXCISION REPAIR`
## NULL
## 
## $`BASE EXCISION REPAIR AP SITE FORMATION`
## NULL
## 
## $`DISEASES OF BASE EXCISION REPAIR`
## NULL
## 
## $`DISEASES OF DNA REPAIR`
## NULL
## 
## $`DISEASES OF MISMATCH REPAIR MMR`
## NULL
## 
## $`DNA DOUBLE STRAND BREAK REPAIR`
## NULL
## 
## $`DNA REPAIR`
## NULL
## 
## $`GAP FILLING DNA REPAIR SYNTHESIS AND LIGATION IN GG NER`
## NULL
## 
## $`GLOBAL GENOME NUCLEOTIDE EXCISION REPAIR GG NER`
## NULL
## 
## $`HOMOLOGY DIRECTED REPAIR`
## NULL
## 
## $`MISMATCH REPAIR`
## NULL
## 
## $`NUCLEOTIDE EXCISION REPAIR`
## NULL
## 
## $`PCNA DEPENDENT LONG PATCH BASE EXCISION REPAIR`
## NULL
## 
## $`POLB DEPENDENT LONG PATCH BASE EXCISION REPAIR`
## NULL
## 
## $`PROTEIN REPAIR`
## NULL
## 
## $`SUMOYLATION OF DNA DAMAGE RESPONSE AND REPAIR PROTEINS`
## NULL
## 
## $`TP53 REGULATES TRANSCRIPTION OF DNA REPAIR GENES`
## NULL
## 
## $`TRANSCRIPTION COUPLED NUCLEOTIDE EXCISION REPAIR TC NER`
## NULL
## 
## $`MITOCHONDRIAL BIOGENESIS`
## NULL
## 
## $`MITOCHONDRIAL CALCIUM ION TRANSPORT`
## NULL
## 
## $`MITOCHONDRIAL FATTY ACID BETA OXIDATION`
## NULL
## 
## $`MITOCHONDRIAL FATTY ACID BETA OXIDATION OF SATURATED FATTY ACIDS`
## NULL
## 
## $`MITOCHONDRIAL FATTY ACID BETA OXIDATION OF UNSATURATED FATTY ACIDS`
## NULL
## 
## $`MITOCHONDRIAL IRON SULFUR CLUSTER BIOGENESIS`
## NULL
## 
## $`MITOCHONDRIAL PROTEIN IMPORT`
## NULL
## 
## $`MITOCHONDRIAL TRANSLATION`
## NULL
## 
## $`MITOCHONDRIAL TRNA AMINOACYLATION`
## NULL
## 
## $`MITOCHONDRIAL UNCOUPLING`
## NULL
## 
## $`TRANSCRIPTIONAL ACTIVATION OF MITOCHONDRIAL BIOGENESIS`
## NULL

Session information

sessionInfo()
## R version 4.1.3 (2022-03-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] ggplot2_3.3.5               beeswarm_0.4.0             
##  [3] kableExtra_1.3.4            eulerr_6.1.1               
##  [5] mitch_1.4.1                 gplots_3.1.1               
##  [7] DESeq2_1.32.0               SummarizedExperiment_1.22.0
##  [9] Biobase_2.52.0              MatrixGenerics_1.4.3       
## [11] matrixStats_0.61.0          GenomicRanges_1.44.0       
## [13] GenomeInfoDb_1.28.4         IRanges_2.26.0             
## [15] S4Vectors_0.30.2            BiocGenerics_0.38.0        
## [17] reshape2_1.4.4             
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.0-3       ellipsis_0.3.2         XVector_0.32.0        
##   [4] rstudioapi_0.13        farver_2.1.0           bit64_4.0.5           
##   [7] AnnotationDbi_1.54.1   fansi_1.0.2            xml2_1.3.3            
##  [10] splines_4.1.3          cachem_1.0.6           geneplotter_1.70.0    
##  [13] knitr_1.37             polyclip_1.10-0        jsonlite_1.8.0        
##  [16] annotate_1.70.0        png_0.1-7              shiny_1.7.1           
##  [19] compiler_4.1.3         httr_1.4.2             assertthat_0.2.1      
##  [22] Matrix_1.4-0           fastmap_1.1.0          cli_3.2.0             
##  [25] later_1.3.0            htmltools_0.5.2        tools_4.1.3           
##  [28] gtable_0.3.0           glue_1.6.2             GenomeInfoDbData_1.2.6
##  [31] dplyr_1.0.8            Rcpp_1.0.8.2           jquerylib_0.1.4       
##  [34] vctrs_0.3.8            Biostrings_2.60.2      svglite_2.1.0         
##  [37] xfun_0.30              polylabelr_0.2.0       stringr_1.4.0         
##  [40] rvest_1.0.2            mime_0.12              lifecycle_1.0.1       
##  [43] gtools_3.9.2           XML_3.99-0.9           zlibbioc_1.38.0       
##  [46] MASS_7.3-55            scales_1.1.1           promises_1.2.0.1      
##  [49] RColorBrewer_1.1-2     yaml_2.3.5             memoise_2.0.1         
##  [52] gridExtra_2.3          sass_0.4.0             reshape_0.8.8         
##  [55] stringi_1.7.6          RSQLite_2.2.10         highr_0.9             
##  [58] genefilter_1.74.1      caTools_1.18.2         BiocParallel_1.26.2   
##  [61] echarts4r_0.4.3        rlang_1.0.2            pkgconfig_2.0.3       
##  [64] systemfonts_1.0.4      bitops_1.0-7           evaluate_0.15         
##  [67] lattice_0.20-45        purrr_0.3.4            labeling_0.4.2        
##  [70] htmlwidgets_1.5.4      bit_4.0.4              tidyselect_1.1.2      
##  [73] GGally_2.1.2           plyr_1.8.6             magrittr_2.0.2        
##  [76] R6_2.5.1               generics_0.1.2         DelayedArray_0.18.0   
##  [79] DBI_1.1.2              withr_2.5.0            pillar_1.7.0          
##  [82] survival_3.3-1         KEGGREST_1.32.0        RCurl_1.98-1.6        
##  [85] tibble_3.1.6           crayon_1.5.0           KernSmooth_2.23-20    
##  [88] utf8_1.2.2             rmarkdown_2.13         locfit_1.5-9.5        
##  [91] grid_4.1.3             blob_1.2.2             digest_0.6.29         
##  [94] webshot_0.5.2          xtable_1.8-4           httpuv_1.6.5          
##  [97] munsell_0.5.0          viridisLite_0.4.0      bslib_0.3.1           
## [100] tcltk_4.1.3