Introduction

suppressPackageStartupMessages({
    library("RhpcBLASctl")
    library("reshape2")
    library("DESeq2")
    library("gplots")
    library("mitch")
    library("eulerr")
    library("limma")
    library("topconfects")
    library("beeswarm")
    library("kableExtra")

    library("dplyr")
    library("network")
})

RhpcBLASctl::blas_set_num_threads(1)

Functions for network chart.

map2color <- function(x,pal,limits=NULL){
  if (is.null(limits)) limits <- range(x)
  pal[findInterval(x,seq(limits[1],limits[2],length.out=length(pal)+1), all.inside=TRUE)]
}

gs2net <- function(gset,em,colfunc=colorRampPalette(c("blue", "white","red"))(n=100)){
  gset <- gset[order(names(gset))]
  mydf <- dplyr::bind_rows(lapply(gset, as.data.frame.list))
  rownames(mydf) <- names(gset)
  j <- apply(mydf,1,function(x) {
    apply(mydf,1,function(y) {
      length(intersect(x,y) ) / length(union(x,y))
    })
  })
  j[lower.tri(j)] <- NA
  j[lower.tri(j,diag=TRUE)] <- 0
  jl <- reshape2::melt(j)
  jl <- jl[which(jl$Var1 != jl$Var2),]
  jl <- jl[which(jl$value != 1),]
  jl <- jl[order(-jl$value),]
  jl <- head(jl,length(gset)*2)
  jl$edgeSize = jl$value/sum(jl$value)
  nodes <- sort(union(jl[,1],jl[,2]))
  lengths <- unlist(lapply(gset,length))
  lengths <- lengths[names(lengths) %in% nodes]
  nodes <- data.frame("nodes"=nodes,"lengths"=lengths)
  nodes$vertexsize <- sqrt(nodes$lengths/sum(nodes$lengths) * 100)
  nodes$es <- em[match(nodes$nodes,em$set),"s.dist"]
  nodes$colours <- map2color(nodes$es,colfunc)
  jl2 <- apply(jl[,1:2],2,as.character)
  jlnet <- network(jl2)
  jlnet$val <- lapply(jlnet$val, function(x) {
    vn <- x[[2]]
    vn <- substr(vn, 1, 60)
    if (nchar(vn) == 60 ) {
      vn <- paste(vn,"...",sep="")
    }
  x[[2]] <- vn
  return(x)
  } )
  plot(jlnet, displaylabels = TRUE, label.col = "steelblue",
       edge.lwd = c(jl$edgeSize) * 100,
       arrowhead.cex = 0,
       label.cex = 0.8, vertex.border = "white",vertex.cex = nodes$vertexsize,
       vertex.col = nodes$colours, edge.col = "black")
  E1 <- min(nodes$es)
  E5 <- max(nodes$es)
  E3 <- mean(c(E1,E5))
  EE <- c(E1,E3,E5)
  legcols <- map2color(EE,colfunc)
  legend("topleft", legend=signif(EE,2) ,title="ES",box.lty=0,
    fill=legcols, cex=0.8)
  S1 <- min(nodes$vertexsize)
  FRAG <- S1/10
  S5 <- max(nodes$vertexsize)
  S3 <- mean(c(S1,S5))
  SS <- c(S1-FRAG,0,S5-FRAG)
  L1 <- min(nodes$lengths)
  L5 <- max(nodes$lengths)
  LL <- paste(" ",c(L1,"",L5))
  legend("topright", legend=LL ,title="no. genes",box.lty=0,
    pt.cex=SS*1, cex=0.9 , pch=19,col="darkgray")
  J1 <- min(jl$edgeSize)
  FRAG <- J1*3
  J5 <- max(jl$edgeSize)
  J3 <- mean(c(J1,J5))
  JJ <- c(J1,J3,J5)
  JL <- JJ+FRAG
  legend("bottomleft", legend=signif(JJ,2) , lwd=JL*50, title="Jaccard",
    box.lty=0, cex=0.9 , lty=1, col="black")
}

networkplot <- function(eres,FDR=0.05,n_sets=20) {
  scores <- eres[[1]][,1]
  names(scores) <- rownames(eres[[1]])
  gs <- eres[[2]]
  eres <- eres[[4]]
  up <- head(eres[eres$p.adjustANOVA < FDR & eres$s.dist > 0,],n_sets)
  n_up <- nrow(up)
  if (n_up >= 5) {
    up_gs <- up[,1]
    up_gs <- gs[which(names(gs) %in% up_gs)]
    topgs_up <- lapply(seq(from=1,to=length(up_gs)),function(i) {
      gsname <- names(up_gs)[i]
      genes <- up_gs[[i]]
      gene_scores <- scores[which(names(scores) %in% genes)]
      top_genes <- names(which(gene_scores > quantile(gene_scores,c(2/3))))
      return(top_genes)
    })
    names(topgs_up) <- names(up_gs)
    gs2net(gset=topgs_up,em=eres,colfunc=colorRampPalette(c("pink","darkred"))(n=100))
  } else {
    message("Can't plot upregulated sets. Fewer than 5 found.")
  }

  dn <- head(eres[eres$p.adjustANOVA < FDR & eres$s.dist < 0,],n_sets)
  n_dn <- nrow(dn)
  if (n_dn >= 5) {
    dn_gs <- dn[,1]
    dn_gs <- gs[which(names(gs) %in% dn_gs)]
    topgs_dn <- lapply(seq(from=1,to=length(dn_gs)),function(i) {
      gsname <- names(dn_gs)[i]
      genes <- dn_gs[[i]]
      gene_scores <- scores[which(names(scores) %in% genes)]
      top_genes <- names(which(gene_scores < quantile(gene_scores,c(1/3))))
      return(top_genes)
    })
    names(topgs_dn) <- names(dn_gs)
    gs2net(gset=topgs_dn,em=eres,colfunc=colorRampPalette(c("darkblue","lightblue"))(n=100))
  } else {
    message("Can't plot downregulated sets. Fewer than 5 found.")
  }
}

network_genes <- function(eres,FDR=0.05,n_sets=20) {
  scores <- eres[[1]][,1]
  names(scores) <- rownames(eres[[1]])
  gs <- eres[[2]]
  eres <- eres[[4]]

  up <- head(eres[eres$p.adjustANOVA < FDR & eres$s.dist > 0,],n_sets)
  n_up <- nrow(up)
  if (n_up > 0) {
    up_gs <- up[,1]
    up_gs <- gs[which(names(gs) %in% up_gs)]
    topgs_up <- lapply(seq(from=1,to=length(up_gs)),function(i) {
      gsname <- names(up_gs)[i]
      genes <- up_gs[[i]]
      gene_scores <- scores[which(names(scores) %in% genes)]
      top_genes <- names(which(gene_scores>2))
      return(top_genes)
    })
    names(topgs_up) <- names(up_gs)
  } else {
    topgs_up <- NULL
    message("No significant upregulated sets to show.")
  }

  dn <- head(eres[eres$p.adjustANOVA < FDR & eres$s.dist < 0,],n_sets)
  n_dn <- nrow(dn)
  if (n_dn > 0) {
    dn_gs <- dn[,1]
    dn_gs <- gs[which(names(gs) %in% dn_gs)]
    topgs_dn <- lapply(seq(from=1,to=length(dn_gs)),function(i) {
      gsname <- names(dn_gs)[i]
      genes <- dn_gs[[i]]
      gene_scores <- scores[which(names(scores) %in% genes)]
      top_genes <- names(which(gene_scores<2))
      return(top_genes)
    })
    names(topgs_dn) <- names(dn_gs)
  } else {
    topgs_dn <- NULL
    message("No significant upregulated sets to show.")
  }
  return(list(c("UP genesets"=topgs_up,"DOWN genesets"=topgs_dn)))
}

Import read counts

Rename samples. Don’t use initials.

tmp <- read.table("raw_data/rna/3col.tsv.gz",header=F)
x <- as.matrix(acast(tmp, V2~V1, value.var="V3", fun.aggregate = sum))
x <- as.data.frame(x)
accession <- sapply((strsplit(rownames(x),"\\|")),"[[",2)
symbol<-sapply((strsplit(rownames(x),"\\|")),"[[",6)
x$geneid <- paste(accession,symbol)

xx <- aggregate(. ~ geneid,x,sum)
rownames(xx) <- xx$geneid
colnames <- gsub("T0R","T0",colnames(xx))
xx$geneid = NULL
xx <- round(xx)
dim(xx)
## [1] 63086    12
colSums(xx)
##  AD-pos_mRNA  AD-pre_mRNA AKH-pos_mRNA AKH-pre_mRNA  AY-pos_mRNA  AY-pre_mRNA 
##     28258815     20506773     22644894     19847636     17583098     11732777 
##  RH-pos_mRNA  RH-pre_mRNA  SB-pos_mRNA  SB-pre_mRNA  ST-pos_mRNA  ST-pre_mRNA 
##     12610349     11542857     14958605      7951981     24373462     20436123
colnames(xx) <- c("pos1","pre1","pos2","pre2","pos3","pre3","pos4","pre4","pos5","pre5","pos6","pre6")

rpm <- apply(xx, 2 , function(x) { x / sum(x) } ) * 1000000
rpm <- rpm[rowMeans(rpm) > 1,]
dim(rpm)
## [1] 14798    12
# gene table
gt <- as.data.frame(rownames(xx))
gt$genesymbol <- sapply(strsplit(gt[,1]," "),"[[",2)
write.table(x=gt,file="gt.tsv",sep="\t")

MDS Plot

mds <- cmdscale(dist(t(xx)))
cols=rep(c("pink","lightblue"),6)
XMIN=min(mds[,1])*1.1
XMAX=max(mds[,1])*1.1
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2", xlim=c(XMIN,XMAX), cex=2,col=cols,pch=19,main="RNA expression", bty="n")
text(mds, labels=colnames(xx) ,cex=1.1)

pdf("rna_mds.pdf",width=5,height=5)
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2", xlim=c(XMIN,XMAX), cex=2,col=cols,pch=19,main="RNA expression", bty="n")
text(mds, labels=colnames(xx) ,cex=1.1)
dev.off()
## png 
##   2

DESeq2

Sample sheet then differential analysis.

Remove genes with less than 10 reads per sample on average.

ss <- as.data.frame(colnames(xx))
ss$pos <- factor(grepl("pos",ss[,1]))
ss$participant <- c("1","1","2","2","3","3","4","4","5","5","6","6")

dim(xx)
## [1] 63086    12
xx <- xx[rowMeans(xx)>=10,]
dim(xx)
## [1] 16971    12
dds <- DESeqDataSetFromMatrix(countData = xx , colData = ss, design = ~ pos )
## converting counts to integer mode
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge_unpaired <- as.data.frame(zz[order(zz$pvalue),])
write.table(dge_unpaired,file="rna_deseq_unpaired.tsv",quote=F,sep="\t")

head(dge_unpaired) |>
  kbl(caption = "Top significant genes in unpaired analysis)") |>
  kable_paper("hover", full_width = F)
Top significant genes in unpaired analysis)
baseMean log2FoldChange lfcSE stat pvalue padj pos1 pre1 pos2 pre2 pos3 pre3 pos4 pre4 pos5 pre5 pos6 pre6
ENSG00000204389.10 HSPA1A 9661.87862 3.220649 0.2248106 14.326056 0 0 14.509585 11.327110 14.313720 11.079875 13.372655 10.824905 14.662918 10.757834 13.943210 11.235981 13.276723 10.746699
ENSG00000204388.7 HSPA1B 7740.52529 3.753745 0.3475392 10.800928 0 0 14.321364 10.730246 13.973664 9.713970 13.051475 10.168243 14.645100 10.065739 13.277423 10.578335 12.788593 10.318778
ENSG00000059804.17 SLC2A3 261.28959 2.981121 0.3029849 9.839175 0 0 9.977001 7.365156 9.070578 7.704917 9.109580 8.012778 9.432930 7.643719 9.172984 7.515668 9.038002 7.949452
ENSG00000132002.9 DNAJB1 2499.31810 1.947194 0.2314796 8.411945 0 0 12.667979 10.418841 12.286233 10.028874 11.334870 10.148086 12.380905 10.422020 11.657551 10.327622 11.200742 10.139620
ENSG00000108551.5 RASD1 76.23498 2.440256 0.3047561 8.007243 0 0 8.193445 7.193235 8.201182 7.517727 8.289186 7.171245 8.345547 7.141452 8.241378 7.515668 7.931817 7.407349
ENSG00000127528.6 KLF2 462.40541 1.654396 0.2102809 7.867552 0 0 9.863643 8.425862 9.585435 8.278253 9.855610 8.669730 9.882036 8.430016 9.838454 8.654425 9.756367 9.165704 
nrow(subset(dge_unpaired,padj<0.05 & log2FoldChange > 0))
## [1] 129
nrow(subset(dge_unpaired,padj<0.05 & log2FoldChange < 0))
## [1] 28
dds <- DESeqDataSetFromMatrix(countData = xx , colData = ss, design = ~ participant + pos )
## converting counts to integer mode
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z <- results(res)
vsd <- vst(dds, blind=FALSE)
zz <- cbind(as.data.frame(z),assay(vsd))
dge <- as.data.frame(zz[order(zz$pvalue),])
write.table(dge,file="rna_deseq.tsv",quote=F,sep="\t")

head(dge) |>
  kbl(caption = "Top significant genes in paired analysis)") |>
  kable_paper("hover", full_width = F)
Top significant genes in paired analysis)
baseMean log2FoldChange lfcSE stat pvalue padj pos1 pre1 pos2 pre2 pos3 pre3 pos4 pre4 pos5 pre5 pos6 pre6
ENSG00000204389.10 HSPA1A 9661.8786 3.143099 0.2161454 14.54160 0 0 14.56451 11.730742 14.376338 11.542009 13.48906 11.353993 14.71246 11.305680 14.023319 11.660457 13.40058 11.297706
ENSG00000265142.9 MIR133A1HG 2355.8815 -1.028605 0.0795560 -12.92933 0 0 11.18380 11.981908 11.604533 12.227695 11.18469 11.756484 11.02525 11.935191 11.225908 11.814978 11.57239 12.390697
ENSG00000119508.18 NR4A3 457.5074 4.072272 0.3594147 11.33029 0 0 11.75432 9.591278 9.796772 9.324578 10.82203 9.528772 10.56108 9.488176 9.973593 9.467974 10.96985 9.587088
ENSG00000127528.6 KLF2 462.4054 1.727463 0.1572805 10.98333 0 0 10.70828 9.918147 10.539975 9.846831 10.70331 10.039391 10.71969 9.920176 10.692706 10.031649 10.64242 10.300906
ENSG00000120694.20 HSPH1 890.8400 1.141477 0.1086403 10.50694 0 0 11.27003 10.500918 10.901201 10.418029 10.98860 10.410880 11.32874 10.914689 11.101868 10.339796 10.94273 10.435839
ENSG00000204388.7 HSPA1B 7740.5253 3.635689 0.3595030 10.11310 0 0 14.38366 11.285948 14.052175 10.616727 13.19459 10.902000 14.69524 10.835690 13.401222 11.178769 12.95765 11.001460 
nrow(subset(dge,padj<0.05 & log2FoldChange > 0))
## [1] 472
nrow(subset(dge,padj<0.05 & log2FoldChange < 0))
## [1] 165

Smear plot

sig <- subset(dge,padj<0.05)
NSIG=nrow(sig)
NDOWN=nrow(subset(sig,log2FoldChange<0))
NUP=nrow(subset(sig,log2FoldChange>0))
NTOT=nrow(dge)
HEADER=paste(NTOT,"genes detected;",NSIG,"@5%FDR;",NUP,"up",NDOWN,"down")
plot(log10(dge$baseMean),dge$log2FoldChange,cex=0.5,pch=19,col="darkgray",
  xlab="log10 basemean",ylab="log2 fold change",main="RNA expression")
points(log10(sig$baseMean),sig$log2FoldChange,cex=0.5,pch=19,col="black")
mtext(HEADER)

pdf("rna_smear.pdf",width=5,height=5)
plot(log10(dge$baseMean),dge$log2FoldChange,cex=0.5,pch=19,col="darkgray",
  xlab="log10 basemean",ylab="log2 fold change",main="RNA expression")
points(log10(sig$baseMean),sig$log2FoldChange,cex=0.5,pch=19,col="black")
mtext(HEADER)
dev.off()
## png 
##   2

Mitch

pw <- gmt_import("ref/c5.go.v2024.1.Hs.symbols.gmt")

names(pw) <- gsub("_", " ", names(pw))

gt <- as.data.frame(rownames(dge))
gt$genesymbol <- sapply(strsplit(gt[,1]," "),"[[",2)
m <- mitch_import(x=dge,DEtype="deseq2",geneTable=gt)
## The input is a single dataframe; one contrast only. Converting
##         it to a list for you.
## Note: Mean no. genes in input = 16971
## Note: no. genes in output = 16923
## Note: estimated proportion of input genes in output = 0.997
head(m)
##                    x
## 5_8S_rRNA  2.8723535
## 7SK        0.1300210
## A1BG       0.2484592
## A2M        3.4477904
## A2M-AS1   -0.5669178
## A4GALT     0.1378113
mres <- mitch_calc(x=m,genesets=pw,minsetsize=5,cores=8,priority="effect")
## Note: Enrichments with large effect sizes may not be
##             statistically significant.
head(subset(mres$enrichment_result,p.adjustANOVA<0.05 & s.dist>0),50) |>
  kbl(caption = "Top upregulated GO terms based on effect size (FDR<0.05)") |>
  kable_paper("hover", full_width = F)
Top upregulated GO terms based on effect size (FDR
set setSize pANOVA s.dist p.adjustANOVA
1145 GOBP ENDOCARDIUM MORPHOGENESIS 5 0.0004658 0.9036529 0.0033947
8446 GOMF NUCLEAR GLUCOCORTICOID RECEPTOR BINDING 11 0.0000005 0.8776125 0.0000076
7757 GOMF C3HC4 TYPE RING FINGER DOMAIN BINDING 5 0.0009659 0.8521811 0.0061947
5126 GOBP REGULATION OF LYMPHANGIOGENESIS 5 0.0012409 0.8338574 0.0076217
5113 GOBP REGULATION OF LIPOPHAGY 5 0.0014736 0.8210900 0.0087528
6323 GOBP TELOMERASE HOLOENZYME COMPLEX ASSEMBLY 6 0.0005975 0.8091663 0.0041868
5284 GOBP REGULATION OF NITRIC OXIDE MEDIATED SIGNAL TRANSDUCTION 5 0.0020478 0.7961461 0.0115809
1967 GOBP LUNG VASCULATURE DEVELOPMENT 5 0.0022720 0.7881310 0.0126210
2618 GOBP NEGATIVE REGULATION OF ENDOTHELIAL CELL DIFFERENTIATION 7 0.0003462 0.7808668 0.0026570
4247 GOBP POSITIVE REGULATION OF WOUND HEALING SPREADING OF EPIDERMAL CELLS 5 0.0025669 0.7786263 0.0139188
5059 GOBP REGULATION OF INNER EAR AUDITORY RECEPTOR CELL DIFFERENTIATION 5 0.0025973 0.7777042 0.0140266
1708 GOBP INHIBITION OF NEUROEPITHELIAL CELL DIFFERENTIATION 6 0.0010019 0.7755315 0.0063894
4396 GOBP PROTEIN LOCALIZATION TO BICELLULAR TIGHT JUNCTION 5 0.0028110 0.7714860 0.0149656
7865 GOMF CO SMAD BINDING 8 0.0001663 0.7686964 0.0014160
3543 GOBP POSITIVE REGULATION OF ARTERY MORPHOGENESIS 6 0.0011556 0.7660145 0.0071835
4665 GOBP REGULATION OF ARTERY MORPHOGENESIS 6 0.0011556 0.7660145 0.0071835
1290 GOBP ESTABLISHMENT OR MAINTENANCE OF CYTOSKELETON POLARITY 5 0.0032951 0.7588604 0.0168899
5601 GOBP REGULATION OF TERMINATION OF DNA TEMPLATED TRANSCRIPTION 5 0.0034472 0.7552429 0.0175212
5698 GOBP REGULATION OF WOUND HEALING SPREADING OF EPIDERMAL CELLS 9 0.0000891 0.7541944 0.0008332
1076 GOBP DORSAL AORTA DEVELOPMENT 9 0.0000996 0.7490048 0.0009164
2952 GOBP NEGATIVE REGULATION OF SPROUTING ANGIOGENESIS 10 0.0000423 0.7476024 0.0004262
3013 GOBP NEGATIVE REGULATION OF T HELPER 17 TYPE IMMUNE RESPONSE 7 0.0006802 0.7414620 0.0046445
4393 GOBP PROTEIN LOCALIZATION TO ADHERENS JUNCTION 6 0.0016881 0.7402416 0.0098012
4253 GOBP POSTSYNAPTIC CYTOSKELETON ORGANIZATION 6 0.0017855 0.7363599 0.0102910
5660 GOBP REGULATION OF T HELPER 17 CELL DIFFERENTIATION 15 0.0000008 0.7359672 0.0000125
2634 GOBP NEGATIVE REGULATION OF ESTABLISHMENT OF PROTEIN LOCALIZATION TO MITOCHONDRION 5 0.0043778 0.7358317 0.0213899
3309 GOBP OVULATION FROM OVARIAN FOLLICLE 7 0.0007491 0.7356856 0.0050167
1077 GOBP DORSAL AORTA MORPHOGENESIS 8 0.0003255 0.7337275 0.0025264
2177 GOBP METANEPHRIC NEPHRON TUBULE MORPHOGENESIS 5 0.0045314 0.7329944 0.0220104
2674 GOBP NEGATIVE REGULATION OF HEMATOPOIETIC STEM CELL DIFFERENTIATION 5 0.0047915 0.7283840 0.0229548
2592 GOBP NEGATIVE REGULATION OF DENDRITE MORPHOGENESIS 5 0.0050050 0.7247665 0.0237893
8146 GOMF HISTONE H4 DEMETHYLASE ACTIVITY 5 0.0051506 0.7223785 0.0243034
1395 GOBP FOREGUT MORPHOGENESIS 7 0.0009540 0.7210249 0.0061356
3978 GOBP POSITIVE REGULATION OF NUCLEAR TRANSCRIBED MRNA CATABOLIC PROCESS DEADENYLATION DEPENDENT DECAY 12 0.0000159 0.7193740 0.0001827
6985 GOCC FILTRATION DIAPHRAGM 5 0.0056835 0.7141270 0.0261529
1191 GOBP ENERGY COUPLED PROTON TRANSMEMBRANE TRANSPORT AGAINST ELECTROCHEMICAL GRADIENT 5 0.0059558 0.7101785 0.0270359
6565 GOBP VENTRICULAR CARDIAC MUSCLE CELL DEVELOPMENT 6 0.0025976 0.7099565 0.0140266
796 GOBP CITRATE METABOLIC PROCESS 5 0.0059992 0.7095638 0.0272000
5119 GOBP REGULATION OF LONG CHAIN FATTY ACID IMPORT ACROSS PLASMA MEMBRANE 6 0.0027990 0.7045970 0.0149103
5548 GOBP REGULATION OF SPONTANEOUS SYNAPTIC TRANSMISSION 5 0.0065665 0.7018797 0.0291961
4911 GOBP REGULATION OF ENDOPLASMIC RETICULUM TUBULAR NETWORK ORGANIZATION 6 0.0029504 0.7007941 0.0155800
7726 GOMF ATP DEPENDENT PROTEIN DISAGGREGASE ACTIVITY 6 0.0029568 0.7006364 0.0156020
3896 GOBP POSITIVE REGULATION OF MICROTUBULE NUCLEATION 8 0.0006246 0.6983447 0.0043396
6561 GOBP VENOUS BLOOD VESSEL MORPHOGENESIS 7 0.0013805 0.6981049 0.0082919
8291 GOMF L LEUCINE BINDING 6 0.0030941 0.6973459 0.0161328
7126 GOCC MICROTUBULE MINUS END 5 0.0070431 0.6958742 0.0307580
3207 GOBP NUCLEOTIDE BINDING OLIGOMERIZATION DOMAIN CONTAINING 1 SIGNALING PATHWAY 7 0.0014464 0.6951660 0.0086208
74 GOBP ADHESION OF SYMBIONT TO HOST CELL 6 0.0032259 0.6943114 0.0166511
8438 GOMF NITRIC OXIDE SYNTHASE REGULATOR ACTIVITY 6 0.0033213 0.6921834 0.0169763
7509 GOCC TRANSCRIPTION FACTOR AP 1 COMPLEX 5 0.0075093 0.6903416 0.0322669 
head(subset(mres$enrichment_result,p.adjustANOVA<0.05 & s.dist<0),50) |>
  kbl(caption = "Top downregulated GO terms based on effect size (FDR<0.05)") |>
  kable_paper("hover", full_width = F)
Top downregulated GO terms based on effect size (FDR
set setSize pANOVA s.dist p.adjustANOVA
6878 GOCC CYTOSOLIC LARGE RIBOSOMAL SUBUNIT 52 0.0000000 -0.8679208 0.0000000
6874 GOCC CYTOPLASMIC SIDE OF ROUGH ENDOPLASMIC RETICULUM MEMBRANE 5 0.0023710 -0.7848209 0.0130590
8240 GOMF LARGE RIBOSOMAL SUBUNIT RRNA BINDING 5 0.0030154 -0.7659298 0.0158131
7766 GOMF CALCIUM DEPENDENT PHOSPHOLIPASE A2 ACTIVITY 5 0.0033264 -0.7581038 0.0169928
6882 GOCC CYTOSOLIC SMALL RIBOSOMAL SUBUNIT 43 0.0000000 -0.7231897 0.0000000
7608 GOMF ACETYLCHOLINE BINDING 5 0.0051258 -0.7227805 0.0242494
7340 GOCC PROTEASOME CORE COMPLEX BETA SUBUNIT COMPLEX 10 0.0001129 -0.7050435 0.0010177
8900 GOMF STRUCTURAL CONSTITUENT OF RIBOSOME 157 0.0000000 -0.6711589 0.0000000
3017 GOBP NEGATIVE REGULATION OF UBIQUITIN PROTEIN LIGASE ACTIVITY 8 0.0010201 -0.6706326 0.0064740
8212 GOMF INTRAMOLECULAR OXIDOREDUCTASE ACTIVITY INTERCONVERTING ALDOSES AND KETOSES 6 0.0044677 -0.6702134 0.0217720
9035 GOMF UBIQUITIN LIGASE INHIBITOR ACTIVITY 9 0.0006626 -0.6553283 0.0045594
6755 GOCC CATSPER COMPLEX 5 0.0114016 -0.6533633 0.0449339
7338 GOCC PROTEASOME CORE COMPLEX 18 0.0000020 -0.6472970 0.0000283
7079 GOCC LARGE RIBOSOMAL SUBUNIT 109 0.0000000 -0.6392663 0.0000000
7609 GOMF ACETYLCHOLINE GATED MONOATOMIC CATION SELECTIVE CHANNEL ACTIVITY 8 0.0025591 -0.6157996 0.0138851
7384 GOCC RIBOSOMAL SUBUNIT 181 0.0000000 -0.6058974 0.0000000
6881 GOCC CYTOSOLIC RIBOSOME 113 0.0000000 -0.6015235 0.0000000
8949 GOMF THREONINE TYPE ENDOPEPTIDASE ACTIVITY 7 0.0075206 -0.5833699 0.0323001
8483 GOMF ODORANT BINDING 8 0.0048302 -0.5753621 0.0230794
7610 GOMF ACETYLCHOLINE RECEPTOR ACTIVITY 8 0.0052371 -0.5700414 0.0245968
7065 GOCC IRON SULFUR CLUSTER ASSEMBLY COMPLEX 12 0.0006726 -0.5668993 0.0046062
6530 GOBP UTP BIOSYNTHETIC PROCESS 8 0.0058770 -0.5623855 0.0267829
7339 GOCC PROTEASOME CORE COMPLEX ALPHA SUBUNIT COMPLEX 7 0.0101835 -0.5608215 0.0413759
4538 GOBP PURINE RIBONUCLEOTIDE SALVAGE 8 0.0065033 -0.5555868 0.0289433
7442 GOCC SMALL RIBOSOMAL SUBUNIT 75 0.0000000 -0.5390298 0.0000000
8428 GOMF NEUROTRANSMITTER RECEPTOR ACTIVITY INVOLVED IN REGULATION OF POSTSYNAPTIC MEMBRANE POTENTIAL 17 0.0001200 -0.5387645 0.0010732
8950 GOMF THREONINE TYPE PEPTIDASE ACTIVITY 12 0.0012457 -0.5381803 0.0076372
8430 GOMF NEUROTRANSMITTER TRANSMEMBRANE TRANSPORTER ACTIVITY 10 0.0040017 -0.5255839 0.0198078
1 GOBP 2FE 2S CLUSTER ASSEMBLY 11 0.0027881 -0.5206631 0.0148611
8664 GOMF POSTSYNAPTIC NEUROTRANSMITTER RECEPTOR ACTIVITY 18 0.0001340 -0.5198988 0.0011778
7143 GOCC MITOCHONDRIAL PROTON TRANSPORTING ATP SYNTHASE COMPLEX COUPLING FACTOR F O 10 0.0050768 -0.5117247 0.0240585
7997 GOMF EXCITATORY EXTRACELLULAR LIGAND GATED MONOATOMIC ION CHANNEL ACTIVITY 17 0.0002678 -0.5105392 0.0021425
7180 GOCC NADPH OXIDASE COMPLEX 8 0.0127535 -0.5084984 0.0491695
1852 GOBP KETONE CATABOLIC PROCESS 11 0.0036348 -0.5063967 0.0183415
4533 GOBP PURINE NUCLEOTIDE SALVAGE 11 0.0038372 -0.5034403 0.0191395
9048 GOMF UBIQUITIN PROTEIN TRANSFERASE INHIBITOR ACTIVITY 12 0.0028508 -0.4973784 0.0150979
1084 GOBP DOUBLE STRAND BREAK REPAIR VIA BREAK INDUCED REPLICATION 9 0.0097980 -0.4971950 0.0401462
968 GOBP DETECTION OF MECHANICAL STIMULUS INVOLVED IN SENSORY PERCEPTION OF SOUND 9 0.0107584 -0.4909542 0.0430333
8484 GOMF OLFACTORY RECEPTOR ACTIVITY 29 0.0000049 -0.4902495 0.0000646
8003 GOMF EXTRACELLULAR LIGAND GATED MONOATOMIC ION CHANNEL ACTIVITY 27 0.0000203 -0.4737479 0.0002271
7657 GOMF ALCOHOL DEHYDROGENASE NADPPLUS ACTIVITY 16 0.0013118 -0.4639942 0.0079586
6028 GOBP RRNA 3 END PROCESSING 10 0.0124583 -0.4563590 0.0482409
7540 GOCC U1 SNRNP 19 0.0006004 -0.4547164 0.0041970
894 GOBP CYTOPLASMIC TRANSLATION 157 0.0000000 -0.4498640 0.0000000
7097 GOCC MANCHETTE 13 0.0058490 -0.4414866 0.0267113
7385 GOCC RIBOSOME 224 0.0000000 -0.4399841 0.0000000
3211 GOBP NUCLEOTIDE SALVAGE 16 0.0025976 -0.4348864 0.0140266
3162 GOBP NUCLEAR MRNA SURVEILLANCE 11 0.0130036 -0.4324847 0.0498425
7352 GOCC PROTON TRANSPORTING ATP SYNTHASE COMPLEX 21 0.0006627 -0.4291574 0.0045594
8427 GOMF NEUROTRANSMITTER RECEPTOR ACTIVITY 25 0.0002228 -0.4265404 0.0018322 
par(mar=c(5,27,3,3))
top <- mres$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)
## [1] 2373
up <- head(subset(top,s.dist>0),20)
dn <- head(subset(top,s.dist<0),20)
top <- rbind(up,dn)
vec=top$s.dist
names(vec)=top$set
names(vec) <- gsub("_"," ",names(vec))
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.65,main="Pre vs Post Exercise",xlab="ES")
grid()

pdf("rna_mitchbar.pdf",width=7,height=7)
par(mar=c(5,27,3,3))
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.65,main="Pre vs Post Exercise",xlab="ES")
grid()
dev.off()
## png 
##   2
par(mar=c(5.1, 4.1, 4.1, 2.1))

if (! file.exists("rna_mitch.html") ) {
  mitch_report(mres,"rna_mitch.html")
}

Network chart and genes.

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

networkplot(mres,FDR=0.05,n_sets=20)

network_genes(mres,FDR=0.05,n_sets=20)
## [[1]]
## [[1]]$`UP genesets.GOBP DORSAL AORTA DEVELOPMENT`
## [1] "DLL4"   "ENG"    "NKX3-1" "RBPJ"   "SRF"   
## 
## [[1]]$`UP genesets.GOBP ENDOCARDIUM MORPHOGENESIS`
## [1] "NOTCH1" "RBPJ"   "SOX17"  "SOX18" 
## 
## [[1]]$`UP genesets.GOBP ESTABLISHMENT OR MAINTENANCE OF CYTOSKELETON POLARITY`
## [1] "AQP1"     "ARHGAP35" "CKAP5"    "LMNA"    
## 
## [[1]]$`UP genesets.GOBP INHIBITION OF NEUROEPITHELIAL CELL DIFFERENTIATION`
## [1] "DLL1"   "HES1"   "JAG1"   "MSX1"   "NOTCH1"
## 
## [[1]]$`UP genesets.GOBP LUNG VASCULATURE DEVELOPMENT`
## [1] "BMP2"   "ERRFI1" "VEGFA" 
## 
## [[1]]$`UP genesets.GOBP NEGATIVE REGULATION OF ENDOTHELIAL CELL DIFFERENTIATION`
## [1] "FOXJ2" "JAG1"  "VEGFA" "ZEB1" 
## 
## [[1]]$`UP genesets.GOBP POSITIVE REGULATION OF ARTERY MORPHOGENESIS`
## [1] "EFNB2"  "NFATC3" "NOTCH1"
## 
## [[1]]$`UP genesets.GOBP POSITIVE REGULATION OF WOUND HEALING SPREADING OF EPIDERMAL CELLS`
## [1] "FERMT2" "MTOR"  
## 
## [[1]]$`UP genesets.GOBP PROTEIN LOCALIZATION TO BICELLULAR TIGHT JUNCTION`
## [1] "ACTG1" "F11R" 
## 
## [[1]]$`UP genesets.GOBP REGULATION OF ARTERY MORPHOGENESIS`
## [1] "EFNB2"  "NFATC3" "NOTCH1"
## 
## [[1]]$`UP genesets.GOBP REGULATION OF INNER EAR AUDITORY RECEPTOR CELL DIFFERENTIATION`
## [1] "DLL1"   "FGF2"   "HES1"   "NOTCH1"
## 
## [[1]]$`UP genesets.GOBP REGULATION OF LIPOPHAGY`
## [1] "ADRB2" "HTT"   "SESN2"
## 
## [[1]]$`UP genesets.GOBP REGULATION OF LYMPHANGIOGENESIS`
## [1] "EPHA2" "FOXC1" "VEGFA"
## 
## [[1]]$`UP genesets.GOBP REGULATION OF NITRIC OXIDE MEDIATED SIGNAL TRANSDUCTION`
## [1] "GUCY1A2" "THBS1"   "VEGFA"  
## 
## [[1]]$`UP genesets.GOBP REGULATION OF TERMINATION OF DNA TEMPLATED TRANSCRIPTION`
## [1] "SCAF4" "SCAF8" "WNK1" 
## 
## [[1]]$`UP genesets.GOBP REGULATION OF WOUND HEALING SPREADING OF EPIDERMAL CELLS`
## [1] "CLASP2" "FERMT2" "MTOR"   "PTEN"  
## 
## [[1]]$`UP genesets.GOBP TELOMERASE HOLOENZYME COMPLEX ASSEMBLY`
## [1] "HSP90AA1" "HSP90AB1" "NAF1"    
## 
## [[1]]$`UP genesets.GOMF C3HC4 TYPE RING FINGER DOMAIN BINDING`
## [1] "DNAJA1" "HSPA1A" "HSPA1B" "HSPA8" 
## 
## [[1]]$`UP genesets.GOMF CO SMAD BINDING`
## [1] "USP9X" "USP9Y"
## 
## [[1]]$`UP genesets.GOMF NUCLEAR GLUCOCORTICOID RECEPTOR BINDING`
## [1] "ETS2"  "FKBP4" "NR4A1" "NR4A2" "NR4A3" "NRIP1" "TACC1" "YWHAH"
## 
## [[1]]$`DOWN genesets.GOBP NEGATIVE REGULATION OF UBIQUITIN PROTEIN LIGASE ACTIVITY`
## [1] "BAG2"   "FBXO5"  "MAD2L1" "MAD2L2" "RPL11"  "RPL23"  "RPL5"   "RPS7"  
## 
## [[1]]$`DOWN genesets.GOCC CATSPER COMPLEX`
## [1] "CATSPER2" "CATSPER3" "CATSPERE" "CATSPERG" "HSPA2"   
## 
## [[1]]$`DOWN genesets.GOCC CYTOPLASMIC SIDE OF ROUGH ENDOPLASMIC RETICULUM MEMBRANE`
## [1] "ALG5"  "RPL27" "RPS26" "RPS28" "RPS29"
## 
## [[1]]$`DOWN genesets.GOCC CYTOSOLIC LARGE RIBOSOMAL SUBUNIT`
##  [1] "RPL10"   "RPL10A"  "RPL11"   "RPL12"   "RPL13"   "RPL13A"  "RPL14"  
##  [8] "RPL15"   "RPL17"   "RPL18"   "RPL18A"  "RPL19"   "RPL21"   "RPL22"  
## [15] "RPL23"   "RPL23A"  "RPL24"   "RPL26"   "RPL26L1" "RPL27"   "RPL27A" 
## [22] "RPL28"   "RPL29"   "RPL3"    "RPL30"   "RPL31"   "RPL32"   "RPL34"  
## [29] "RPL35"   "RPL35A"  "RPL36"   "RPL36A"  "RPL36AL" "RPL37"   "RPL37A" 
## [36] "RPL38"   "RPL39"   "RPL3L"   "RPL4"    "RPL41"   "RPL5"    "RPL6"   
## [43] "RPL7"    "RPL7A"   "RPL7L1"  "RPL8"    "RPL9"    "RPLP0"   "RPLP1"  
## [50] "RPLP2"   "UBA52"   "ZCCHC17"
## 
## [[1]]$`DOWN genesets.GOCC CYTOSOLIC RIBOSOME`
##   [1] "ABCE1"   "APOD"    "ASCC2"   "DHX29"   "EEF1A1"  "EIF2A"   "EIF2AK4"
##   [8] "EIF2D"   "ETF1"    "FAU"     "GCN1"    "HBS1L"   "LARP4"   "MCTS1"  
##  [15] "METAP1"  "NEMF"    "PELO"    "RACK1"   "RNF10"   "RNF14"   "RNF25"  
##  [22] "RPL10"   "RPL10A"  "RPL11"   "RPL12"   "RPL13"   "RPL13A"  "RPL14"  
##  [29] "RPL15"   "RPL17"   "RPL18"   "RPL18A"  "RPL19"   "RPL21"   "RPL22"  
##  [36] "RPL23"   "RPL23A"  "RPL24"   "RPL26"   "RPL26L1" "RPL27"   "RPL27A" 
##  [43] "RPL28"   "RPL29"   "RPL3"    "RPL30"   "RPL31"   "RPL32"   "RPL34"  
##  [50] "RPL35"   "RPL35A"  "RPL36"   "RPL36A"  "RPL36AL" "RPL37"   "RPL37A" 
##  [57] "RPL38"   "RPL39"   "RPL3L"   "RPL4"    "RPL41"   "RPL5"    "RPL6"   
##  [64] "RPL7"    "RPL7A"   "RPL7L1"  "RPL8"    "RPL9"    "RPLP0"   "RPLP1"  
##  [71] "RPLP2"   "RPS10"   "RPS11"   "RPS12"   "RPS13"   "RPS14"   "RPS15"  
##  [78] "RPS15A"  "RPS16"   "RPS17"   "RPS18"   "RPS19"   "RPS2"    "RPS20"  
##  [85] "RPS21"   "RPS23"   "RPS24"   "RPS25"   "RPS26"   "RPS27"   "RPS27A" 
##  [92] "RPS27L"  "RPS28"   "RPS29"   "RPS3"    "RPS3A"   "RPS4X"   "RPS4Y1" 
##  [99] "RPS5"    "RPS6"    "RPS7"    "RPS8"    "RPS9"    "RPSA"    "UBA52"  
## [106] "USP10"   "ZCCHC17" "ZNF598" 
## 
## [[1]]$`DOWN genesets.GOCC CYTOSOLIC SMALL RIBOSOMAL SUBUNIT`
##  [1] "DHX29"  "EIF2A"  "EIF2D"  "FAU"    "LARP4"  "MCTS1"  "RACK1"  "RPS10" 
##  [9] "RPS11"  "RPS12"  "RPS13"  "RPS14"  "RPS15"  "RPS15A" "RPS16"  "RPS17" 
## [17] "RPS18"  "RPS19"  "RPS2"   "RPS20"  "RPS21"  "RPS23"  "RPS24"  "RPS25" 
## [25] "RPS26"  "RPS27"  "RPS27A" "RPS27L" "RPS28"  "RPS29"  "RPS3"   "RPS3A" 
## [33] "RPS4X"  "RPS4Y1" "RPS5"   "RPS6"   "RPS7"   "RPS8"   "RPS9"   "RPSA"  
## [41] "UBA52" 
## 
## [[1]]$`DOWN genesets.GOCC LARGE RIBOSOMAL SUBUNIT`
##   [1] "GADD45GIP1" "MPV17L2"    "MRPL1"      "MRPL10"     "MRPL11"    
##   [6] "MRPL12"     "MRPL13"     "MRPL14"     "MRPL15"     "MRPL16"    
##  [11] "MRPL17"     "MRPL18"     "MRPL19"     "MRPL2"      "MRPL20"    
##  [16] "MRPL21"     "MRPL22"     "MRPL23"     "MRPL24"     "MRPL27"    
##  [21] "MRPL28"     "MRPL3"      "MRPL30"     "MRPL32"     "MRPL33"    
##  [26] "MRPL34"     "MRPL35"     "MRPL36"     "MRPL37"     "MRPL38"    
##  [31] "MRPL4"      "MRPL40"     "MRPL41"     "MRPL42"     "MRPL43"    
##  [36] "MRPL44"     "MRPL45"     "MRPL46"     "MRPL47"     "MRPL48"    
##  [41] "MRPL49"     "MRPL50"     "MRPL51"     "MRPL52"     "MRPL53"    
##  [46] "MRPL54"     "MRPL55"     "MRPL57"     "MRPL58"     "MRPL9"     
##  [51] "MRPS18A"    "MRPS30"     "MTERF4"     "NSUN3"      "NSUN4"     
##  [56] "RBM3"       "RPL10"      "RPL10A"     "RPL11"      "RPL12"     
##  [61] "RPL13"      "RPL13A"     "RPL14"      "RPL15"      "RPL17"     
##  [66] "RPL18"      "RPL18A"     "RPL19"      "RPL21"      "RPL22"     
##  [71] "RPL23"      "RPL23A"     "RPL24"      "RPL26"      "RPL26L1"   
##  [76] "RPL27"      "RPL27A"     "RPL28"      "RPL29"      "RPL3"      
##  [81] "RPL30"      "RPL31"      "RPL32"      "RPL34"      "RPL35"     
##  [86] "RPL35A"     "RPL36"      "RPL36A"     "RPL36AL"    "RPL37"     
##  [91] "RPL37A"     "RPL38"      "RPL39"      "RPL3L"      "RPL4"      
##  [96] "RPL41"      "RPL5"       "RPL6"       "RPL7"       "RPL7A"     
## [101] "RPL7L1"     "RPL8"       "RPL9"       "RPLP0"      "RPLP1"     
## [106] "RPLP2"      "UBA52"      "ZCCHC17"   
## 
## [[1]]$`DOWN genesets.GOCC PROTEASOME CORE COMPLEX`
##  [1] "PSMA1"  "PSMA2"  "PSMA3"  "PSMA4"  "PSMA5"  "PSMA6"  "PSMA7"  "PSMB1" 
##  [9] "PSMB10" "PSMB2"  "PSMB3"  "PSMB4"  "PSMB5"  "PSMB6"  "PSMB7"  "PSMB8" 
## [17] "PSMB9"  "PSMF1" 
## 
## [[1]]$`DOWN genesets.GOCC PROTEASOME CORE COMPLEX BETA SUBUNIT COMPLEX`
##  [1] "PSMB1"  "PSMB10" "PSMB2"  "PSMB3"  "PSMB4"  "PSMB5"  "PSMB6"  "PSMB7" 
##  [9] "PSMB8"  "PSMB9" 
## 
## [[1]]$`DOWN genesets.GOCC RIBOSOMAL SUBUNIT`
##   [1] "AURKAIP1"   "CHCHD1"     "DAP3"       "DHX29"      "EIF2A"     
##   [6] "EIF2D"      "FAU"        "GADD45GIP1" "LARP4"      "MCTS1"     
##  [11] "MPV17L2"    "MRPL1"      "MRPL10"     "MRPL11"     "MRPL12"    
##  [16] "MRPL13"     "MRPL14"     "MRPL15"     "MRPL16"     "MRPL17"    
##  [21] "MRPL18"     "MRPL19"     "MRPL2"      "MRPL20"     "MRPL21"    
##  [26] "MRPL22"     "MRPL23"     "MRPL24"     "MRPL27"     "MRPL28"    
##  [31] "MRPL3"      "MRPL30"     "MRPL32"     "MRPL33"     "MRPL34"    
##  [36] "MRPL35"     "MRPL36"     "MRPL37"     "MRPL38"     "MRPL4"     
##  [41] "MRPL40"     "MRPL41"     "MRPL42"     "MRPL43"     "MRPL44"    
##  [46] "MRPL45"     "MRPL46"     "MRPL47"     "MRPL48"     "MRPL49"    
##  [51] "MRPL50"     "MRPL51"     "MRPL52"     "MRPL53"     "MRPL54"    
##  [56] "MRPL55"     "MRPL57"     "MRPL58"     "MRPL9"      "MRPS10"    
##  [61] "MRPS11"     "MRPS12"     "MRPS14"     "MRPS15"     "MRPS16"    
##  [66] "MRPS17"     "MRPS18A"    "MRPS18B"    "MRPS18C"    "MRPS2"     
##  [71] "MRPS21"     "MRPS22"     "MRPS23"     "MRPS24"     "MRPS25"    
##  [76] "MRPS26"     "MRPS27"     "MRPS28"     "MRPS30"     "MRPS31"    
##  [81] "MRPS33"     "MRPS34"     "MRPS35"     "MRPS5"      "MRPS6"     
##  [86] "MRPS7"      "MRPS9"      "MTERF4"     "NSUN3"      "NSUN4"     
##  [91] "PTCD3"      "RACK1"      "RBM3"       "RPL10"      "RPL10A"    
##  [96] "RPL11"      "RPL12"      "RPL13"      "RPL13A"     "RPL14"     
## [101] "RPL15"      "RPL17"      "RPL18"      "RPL18A"     "RPL19"     
## [106] "RPL21"      "RPL22"      "RPL23"      "RPL23A"     "RPL24"     
## [111] "RPL26"      "RPL26L1"    "RPL27"      "RPL27A"     "RPL28"     
## [116] "RPL29"      "RPL3"       "RPL30"      "RPL31"      "RPL32"     
## [121] "RPL34"      "RPL35"      "RPL35A"     "RPL36"      "RPL36A"    
## [126] "RPL36AL"    "RPL37"      "RPL37A"     "RPL38"      "RPL39"     
## [131] "RPL3L"      "RPL4"       "RPL41"      "RPL5"       "RPL6"      
## [136] "RPL7"       "RPL7A"      "RPL7L1"     "RPL8"       "RPL9"      
## [141] "RPLP0"      "RPLP1"      "RPLP2"      "RPS10"      "RPS11"     
## [146] "RPS12"      "RPS13"      "RPS14"      "RPS15"      "RPS15A"    
## [151] "RPS16"      "RPS17"      "RPS18"      "RPS19"      "RPS2"      
## [156] "RPS20"      "RPS21"      "RPS23"      "RPS24"      "RPS25"     
## [161] "RPS26"      "RPS27"      "RPS27A"     "RPS27L"     "RPS28"     
## [166] "RPS29"      "RPS3"       "RPS3A"      "RPS4X"      "RPS4Y1"    
## [171] "RPS5"       "RPS6"       "RPS7"       "RPS8"       "RPS9"      
## [176] "RPSA"       "UBA52"      "ZCCHC17"   
## 
## [[1]]$`DOWN genesets.GOMF ACETYLCHOLINE BINDING`
## [1] "ACHE"   "CHRNA1" "CHRNA7" "CHRNB1" "CHRND" 
## 
## [[1]]$`DOWN genesets.GOMF ACETYLCHOLINE GATED MONOATOMIC CATION SELECTIVE CHANNEL ACTIVITY`
## [1] "CHRNA1"  "CHRNA10" "CHRNA5"  "CHRNA7"  "CHRNB1"  "CHRND"   "CHRNE"  
## [8] "CHRNG"  
## 
## [[1]]$`DOWN genesets.GOMF ACETYLCHOLINE RECEPTOR ACTIVITY`
## [1] "ANXA9"  "CHRNA1" "CHRNA5" "CHRNA7" "CHRNB1" "CHRND"  "CHRNE"  "CHRNG" 
## 
## [[1]]$`DOWN genesets.GOMF CALCIUM DEPENDENT PHOSPHOLIPASE A2 ACTIVITY`
## [1] "PLA2G12A" "PLA2G2C"  "PLA2G4B"  "PLA2G4C"  "PLA2G4F" 
## 
## [[1]]$`DOWN genesets.GOMF INTRAMOLECULAR OXIDOREDUCTASE ACTIVITY INTERCONVERTING ALDOSES AND KETOSES`
## [1] "GPI"  "HYI"  "MPI"  "MRI1" "RPIA" "TPI1"
## 
## [[1]]$`DOWN genesets.GOMF LARGE RIBOSOMAL SUBUNIT RRNA BINDING`
## [1] "MRPL11" "RPL12"  "RPL23"  "RPLP0"  "RPUSD4"
## 
## [[1]]$`DOWN genesets.GOMF ODORANT BINDING`
## [1] "OR10K1" "OR14J1" "OR5A1"  "OR5A2"  "OR5AS1" "OR5AU1" "OR5H14" "OR8A1" 
## 
## [[1]]$`DOWN genesets.GOMF STRUCTURAL CONSTITUENT OF RIBOSOME`
##   [1] "DAP3"      "FAU"       "LINC01004" "MRPL1"     "MRPL10"    "MRPL11"   
##   [7] "MRPL12"    "MRPL13"    "MRPL14"    "MRPL15"    "MRPL16"    "MRPL17"   
##  [13] "MRPL18"    "MRPL19"    "MRPL2"     "MRPL20"    "MRPL21"    "MRPL22"   
##  [19] "MRPL23"    "MRPL24"    "MRPL27"    "MRPL28"    "MRPL3"     "MRPL30"   
##  [25] "MRPL32"    "MRPL33"    "MRPL34"    "MRPL35"    "MRPL36"    "MRPL37"   
##  [31] "MRPL4"     "MRPL41"    "MRPL42"    "MRPL43"    "MRPL45"    "MRPL46"   
##  [37] "MRPL47"    "MRPL49"    "MRPL51"    "MRPL52"    "MRPL54"    "MRPL55"   
##  [43] "MRPL57"    "MRPL9"     "MRPS11"    "MRPS12"    "MRPS14"    "MRPS15"   
##  [49] "MRPS16"    "MRPS17"    "MRPS18A"   "MRPS18B"   "MRPS18C"   "MRPS2"    
##  [55] "MRPS21"    "MRPS22"    "MRPS23"    "MRPS24"    "MRPS25"    "MRPS30"   
##  [61] "MRPS31"    "MRPS33"    "MRPS34"    "MRPS35"    "MRPS5"     "MRPS6"    
##  [67] "MRPS7"     "MRPS9"     "NDUFA7"    "RPL10"     "RPL10A"    "RPL11"    
##  [73] "RPL12"     "RPL13"     "RPL13A"    "RPL14"     "RPL15"     "RPL17"    
##  [79] "RPL18"     "RPL18A"    "RPL19"     "RPL21"     "RPL22"     "RPL22L1"  
##  [85] "RPL23"     "RPL23A"    "RPL24"     "RPL26"     "RPL26L1"   "RPL27"    
##  [91] "RPL27A"    "RPL28"     "RPL29"     "RPL3"      "RPL30"     "RPL31"    
##  [97] "RPL32"     "RPL34"     "RPL35"     "RPL35A"    "RPL36"     "RPL36A"   
## [103] "RPL36AL"   "RPL37"     "RPL37A"    "RPL38"     "RPL39"     "RPL3L"    
## [109] "RPL4"      "RPL41"     "RPL5"      "RPL6"      "RPL7"      "RPL7A"    
## [115] "RPL7L1"    "RPL8"      "RPL9"      "RPLP0"     "RPLP1"     "RPLP2"    
## [121] "RPS10"     "RPS11"     "RPS12"     "RPS13"     "RPS14"     "RPS15"    
## [127] "RPS15A"    "RPS16"     "RPS17"     "RPS18"     "RPS19"     "RPS2"     
## [133] "RPS20"     "RPS21"     "RPS23"     "RPS24"     "RPS25"     "RPS26"    
## [139] "RPS27"     "RPS27A"    "RPS27L"    "RPS28"     "RPS29"     "RPS3"     
## [145] "RPS3A"     "RPS4X"     "RPS4Y1"    "RPS5"      "RPS6"      "RPS7"     
## [151] "RPS8"      "RPS9"      "RPSA"      "RSL24D1"   "SRBD1"     "UBA52"    
## 
## [[1]]$`DOWN genesets.GOMF THREONINE TYPE ENDOPEPTIDASE ACTIVITY`
## [1] "PSMB10" "PSMB5"  "PSMB6"  "PSMB7"  "PSMB8"  "PSMB9"  "TASP1" 
## 
## [[1]]$`DOWN genesets.GOMF UBIQUITIN LIGASE INHIBITOR ACTIVITY`
## [1] "FBXO5" "HTRA2" "RPL11" "RPL23" "RPL37" "RPL5"  "RPS15" "RPS20" "RPS7"

SLC

SLC2A4 <- rpm[grep("ENSG00000181856",rownames(rpm)),]
SLC2A4
##     pos1     pre1     pos2     pre2     pos3     pre3     pos4     pre4 
## 192.7894 211.7837 268.4932 290.2613 230.0505 250.5801 176.6010 200.4703 
##     pos5     pre5     pos6     pre6 
## 180.0970 201.4592 192.6686 252.2494
SLC2A4 <- matrix(SLC2A4,nrow=2)
colnames(SLC2A4) <- 1:6
rownames(SLC2A4) <- c("Post","Pre")
SLC2A4 <- SLC2A4[c(2,1),]
barplot(SLC2A4,beside=TRUE,col=c("#0047AB", "#D2042D"),ylim=c(0,300),
  main="SLC2A4 mRNA",ylab="counts per million",xlab="participant")
legend(11,300, legend=c("Pre", "Post"), fill=c("#0047AB","#D2042D"), cex=1)

pdf("rna_SLC2A4.pdf",width=5,height=5)
barplot(SLC2A4,beside=TRUE,col=c("#0047AB", "#D2042D"),ylim=c(0,300),
  main="SLC2A4 mRNA",ylab="counts per million",xlab="participant")
legend(11,300, legend=c("Pre", "Post"), fill=c("#0047AB","#D2042D"), cex=1)
dev.off()
## png 
##   2
tpm <- apply(x[,1:12], 2 , function(x) { x / sum(x) } ) * 1000000
tpm <- tpm[rowMeans(tpm) > 1,]
dim(tpm)
## [1] 36942    12
SLC2A4 <- tpm[grep("ENSG00000181856",rownames(tpm)),]
rownames(SLC2A4) <- sapply(strsplit(rownames(SLC2A4),"\\|"),"[[",1)
SLC2A4
##                    AD-pos_mRNA AD-pre_mRNA AKH-pos_mRNA AKH-pre_mRNA
## ENST00000317370.13 145.2675550   153.17734   175.516303    203.11333
## ENST00000424875.2    0.6459906     0.00000     1.173132      0.00000
## ENST00000572485.5   45.8908386    57.33124    91.091064     86.04494
##                    AY-pos_mRNA AY-pre_mRNA RH-pos_mRNA RH-pre_mRNA SB-pos_mRNA
## ENST00000317370.13  161.195089  189.143942  133.162401   116.59716  132.028166
## ENST00000424875.2     3.475109    1.603876    1.701429     0.00000    2.563606
## ENST00000572485.5    65.349422   59.305853   41.765433    83.17798   45.505262
##                    SB-pre_mRNA ST-pos_mRNA ST-pre_mRNA
## ENST00000317370.13   158.05147  136.457650   179.10165
## ENST00000424875.2      1.30145    1.856872     2.97331
## ENST00000572485.5     42.03799   54.349366    70.16995

HATs and HDACs

hats <- c("ENSG00000010282.13 HAT1",
  "ENSG00000122390.19 NAA60",
  "ENSG00000108773.11 KAT2A",
  "ENSG000001141663.8 KAT2B",
  "ENSG00000172977.13 KAT5",
  "ENSG00000083168.11 KAT6A",
  "ENSG00000156650.14 KAT6B",
  "ENSG00000136504.15 KAT7",
  "ENSG00000103510.20 KAT8",
  "ENSG00000100393.16 EP300",
  "ENSG00000005339.15 CREBBP",
  "ENSG00000147133.18 TAF1",
  "ENSG00000125484.12 GTF3C4",
  "ENSG00000084676.16 NCOA1",
  "ENSG00000124151.19 NCOA3",
  "ENSG00000140396.13 NCOA2",
  "ENSG00000134852.15 CLOCK")

hdacs <- c("ENSG00000116478.12 HDAC1",
  "ENSG00000196591.12 HDAC2",
  "ENSG00000171720.10 HDAC3",
  "ENSG00000068024.18 HDAC4",
  "ENSG00000108840.17 HDAC5",
  "ENSG00000094631.22 HDAC6",
  "ENSG00000061273.18 HDAC7",
  "ENSG00000147099.21 HDAC8",
  "ENSG00000048052.25 HDAC9",
  "ENSG00000100429.18 HDAC10",
  "ENSG00000163517.15 HDAC11",
  "ENSG00000096717.12 SIRT1",
  "ENSG00000068903.21 SIRT2",
  "ENSG00000142082.15 SIRT3",
  "ENSG00000089163.4 SIRT4",
  "ENSG00000124523.17 SIRT5",
  "ENSG00000077463.15 SIRT6",
  "ENSG00000187531.14 SIRT7")

y1 <- rpm[which(rownames(rpm) %in% hats),,drop=FALSE]

gbarplot <- function(i,y) {
  y <- y[i,,drop=FALSE]
  gname=rownames(y)[1]
  yy <- matrix(y[1,],nrow=2)
  colnames(yy) <- 1:6
  rownames(yy) <- c("Post","Pre")
  yy <- yy[c(2,1),]

  barplot(yy,beside=TRUE,col=c("#0047AB", "#D2042D"),
    main=gname,ylab="counts per million",xlab="participant")
  legend(11,300, legend=c("Pre", "Post"), fill=c("#0047AB","#D2042D"), cex=1)

  legend(11,300, legend=c("Pre", "Post"), fill=c("#0047AB","#D2042D"), cex=1)
}

lapply(1:nrow(y1), gbarplot, y1)

## [[1]]
## [[1]]$rect
## [[1]]$rect$w
## [1] 4.044734
## 
## [[1]]$rect$h
## [1] 14.77587
## 
## [[1]]$rect$left
## [1] 11
## 
## [[1]]$rect$top
## [1] 300
## 
## 
## [[1]]$text
## [[1]]$text$x
## [1] 13.05085 13.05085
## 
## [[1]]$text$y
## [1] 295.0747 290.1494
## 
## 
## 
## [[2]]
## [[2]]$rect
## [[2]]$rect$w
## [1] 4.044734
## 
## [[2]]$rect$h
## [1] 7.293703
## 
## [[2]]$rect$left
## [1] 11
## 
## [[2]]$rect$top
## [1] 300
## 
## 
## [[2]]$text
## [[2]]$text$x
## [1] 13.05085 13.05085
## 
## [[2]]$text$y
## [1] 297.5688 295.1375
## 
## 
## 
## [[3]]
## [[3]]$rect
## [[3]]$rect$w
## [1] 4.044734
## 
## [[3]]$rect$h
## [1] 12.81928
## 
## [[3]]$rect$left
## [1] 11
## 
## [[3]]$rect$top
## [1] 300
## 
## 
## [[3]]$text
## [[3]]$text$x
## [1] 13.05085 13.05085
## 
## [[3]]$text$y
## [1] 295.7269 291.4538
## 
## 
## 
## [[4]]
## [[4]]$rect
## [[4]]$rect$w
## [1] 4.044734
## 
## [[4]]$rect$h
## [1] 10.11755
## 
## [[4]]$rect$left
## [1] 11
## 
## [[4]]$rect$top
## [1] 300
## 
## 
## [[4]]$text
## [[4]]$text$x
## [1] 13.05085 13.05085
## 
## [[4]]$text$y
## [1] 296.6275 293.2550
## 
## 
## 
## [[5]]
## [[5]]$rect
## [[5]]$rect$w
## [1] 4.044734
## 
## [[5]]$rect$h
## [1] 9.072821
## 
## [[5]]$rect$left
## [1] 11
## 
## [[5]]$rect$top
## [1] 300
## 
## 
## [[5]]$text
## [[5]]$text$x
## [1] 13.05085 13.05085
## 
## [[5]]$text$y
## [1] 296.9757 293.9515
## 
## 
## 
## [[6]]
## [[6]]$rect
## [[6]]$rect$w
## [1] 4.044734
## 
## [[6]]$rect$h
## [1] 14.48734
## 
## [[6]]$rect$left
## [1] 11
## 
## [[6]]$rect$top
## [1] 300
## 
## 
## [[6]]$text
## [[6]]$text$x
## [1] 13.05085 13.05085
## 
## [[6]]$text$y
## [1] 295.1709 290.3418
## 
## 
## 
## [[7]]
## [[7]]$rect
## [[7]]$rect$w
## [1] 4.044734
## 
## [[7]]$rect$h
## [1] 7.198351
## 
## [[7]]$rect$left
## [1] 11
## 
## [[7]]$rect$top
## [1] 300
## 
## 
## [[7]]$text
## [[7]]$text$x
## [1] 13.05085 13.05085
## 
## [[7]]$text$y
## [1] 297.6005 295.2011
## 
## 
## 
## [[8]]
## [[8]]$rect
## [[8]]$rect$w
## [1] 4.044734
## 
## [[8]]$rect$h
## [1] 12.77775
## 
## [[8]]$rect$left
## [1] 11
## 
## [[8]]$rect$top
## [1] 300
## 
## 
## [[8]]$text
## [[8]]$text$x
## [1] 13.05085 13.05085
## 
## [[8]]$text$y
## [1] 295.7408 291.4815
## 
## 
## 
## [[9]]
## [[9]]$rect
## [[9]]$rect$w
## [1] 4.044734
## 
## [[9]]$rect$h
## [1] 3.448627
## 
## [[9]]$rect$left
## [1] 11
## 
## [[9]]$rect$top
## [1] 300
## 
## 
## [[9]]$text
## [[9]]$text$x
## [1] 13.05085 13.05085
## 
## [[9]]$text$y
## [1] 298.8505 297.7009
## 
## 
## 
## [[10]]
## [[10]]$rect
## [[10]]$rect$w
## [1] 4.044734
## 
## [[10]]$rect$h
## [1] 6.163766
## 
## [[10]]$rect$left
## [1] 11
## 
## [[10]]$rect$top
## [1] 300
## 
## 
## [[10]]$text
## [[10]]$text$x
## [1] 13.05085 13.05085
## 
## [[10]]$text$y
## [1] 297.9454 295.8908
## 
## 
## 
## [[11]]
## [[11]]$rect
## [[11]]$rect$w
## [1] 4.044734
## 
## [[11]]$rect$h
## [1] 17.55452
## 
## [[11]]$rect$left
## [1] 11
## 
## [[11]]$rect$top
## [1] 300
## 
## 
## [[11]]$text
## [[11]]$text$x
## [1] 13.05085 13.05085
## 
## [[11]]$text$y
## [1] 294.1485 288.2970
## 
## 
## 
## [[12]]
## [[12]]$rect
## [[12]]$rect$w
## [1] 4.044734
## 
## [[12]]$rect$h
## [1] 4.338335
## 
## [[12]]$rect$left
## [1] 11
## 
## [[12]]$rect$top
## [1] 300
## 
## 
## [[12]]$text
## [[12]]$text$x
## [1] 13.05085 13.05085
## 
## [[12]]$text$y
## [1] 298.5539 297.1078
## 
## 
## 
## [[13]]
## [[13]]$rect
## [[13]]$rect$w
## [1] 4.044734
## 
## [[13]]$rect$h
## [1] 5.165018
## 
## [[13]]$rect$left
## [1] 11
## 
## [[13]]$rect$top
## [1] 300
## 
## 
## [[13]]$text
## [[13]]$text$x
## [1] 13.05085 13.05085
## 
## [[13]]$text$y
## [1] 298.2783 296.5567
## 
## 
## 
## [[14]]
## [[14]]$rect
## [[14]]$rect$w
## [1] 4.044734
## 
## [[14]]$rect$h
## [1] 10.52618
## 
## [[14]]$rect$left
## [1] 11
## 
## [[14]]$rect$top
## [1] 300
## 
## 
## [[14]]$text
## [[14]]$text$x
## [1] 13.05085 13.05085
## 
## [[14]]$text$y
## [1] 296.4913 292.9825
## 
## 
## 
## [[15]]
## [[15]]$rect
## [[15]]$rect$w
## [1] 4.044734
## 
## [[15]]$rect$h
## [1] 3.294816
## 
## [[15]]$rect$left
## [1] 11
## 
## [[15]]$rect$top
## [1] 300
## 
## 
## [[15]]$text
## [[15]]$text$x
## [1] 13.05085 13.05085
## 
## [[15]]$text$y
## [1] 298.9017 297.8035
dge1 <- dge[which(rownames(dge) %in% hats),]

head(dge1) |>
  kbl(caption = "HAT genes)") |>
  kable_paper("hover", full_width = F)
HAT genes)
baseMean log2FoldChange lfcSE stat pvalue padj pos1 pre1 pos2 pre2 pos3 pre3 pos4 pre4 pos5 pre5 pos6 pre6
ENSG00000108773.11 KAT2A 823.7213 -0.3904717 0.0927294 -4.210873 0.0000254 0.0014766 10.66647 10.86577 10.55692 10.73603 10.66988 10.91499 10.90893 11.03982 10.64062 11.01623 10.59578 10.68415
ENSG00000124151.19 NCOA3 864.1477 0.2705133 0.0923327 2.929767 0.0033922 0.0560482 10.83693 10.67429 10.79039 10.76885 10.94746 10.69712 10.68816 10.67039 10.99910 10.83093 11.06034 10.81702
ENSG00000005339.15 CREBBP 1002.6364 0.2425909 0.0900621 2.693595 0.0070686 0.0916149 11.04177 10.89736 10.93579 10.90716 11.16227 11.08192 10.90787 10.87034 10.97749 10.68153 10.96481 10.74819
ENSG00000100393.16 EP300 747.8770 0.2721092 0.1034069 2.631441 0.0085024 0.1005081 10.85689 10.67519 10.57045 10.58967 10.86418 10.77953 10.79672 10.70386 10.80649 10.46075 10.76856 10.62672
ENSG00000083168.11 KAT6A 562.2708 0.2586337 0.1076801 2.401870 0.0163115 0.1471017 10.58150 10.40080 10.56966 10.47150 10.57632 10.42757 10.47929 10.51807 10.60171 10.54966 10.63057 10.39218
ENSG00000084676.16 NCOA1 868.3112 0.2277903 0.0969115 2.350499 0.0187482 0.1609802 11.12050 11.01614 10.75664 10.70299 10.88512 10.72777 10.78853 10.80932 10.84696 10.52308 10.86125 10.75817 
HEADER=paste(NTOT,"genes detected;",NSIG,"@5%FDR;",NUP,"up",NDOWN,"down")

plot(log10(dge$baseMean),dge$log2FoldChange,cex=0.5,pch=19,col="darkgray",
  xlab="log10 basemean",ylab="log2 fold change",main="RNA expression")
points(log10(dge1$baseMean),dge1$log2FoldChange,cex=0.5,pch=19,col="black")
mtext("HATs")

y2 <- rpm[which(rownames(rpm) %in% hdacs),,drop=FALSE]

lapply(1:nrow(y2), gbarplot, y2)

## [[1]]
## [[1]]$rect
## [[1]]$rect$w
## [1] 4.044734
## 
## [[1]]$rect$h
## [1] 6.5625
## 
## [[1]]$rect$left
## [1] 11
## 
## [[1]]$rect$top
## [1] 300
## 
## 
## [[1]]$text
## [[1]]$text$x
## [1] 13.05085 13.05085
## 
## [[1]]$text$y
## [1] 297.8125 295.6250
## 
## 
## 
## [[2]]
## [[2]]$rect
## [[2]]$rect$w
## [1] 4.044734
## 
## [[2]]$rect$h
## [1] 4.292697
## 
## [[2]]$rect$left
## [1] 11
## 
## [[2]]$rect$top
## [1] 300
## 
## 
## [[2]]$text
## [[2]]$text$x
## [1] 13.05085 13.05085
## 
## [[2]]$text$y
## [1] 298.5691 297.1382
## 
## 
## 
## [[3]]
## [[3]]$rect
## [[3]]$rect$w
## [1] 4.044734
## 
## [[3]]$rect$h
## [1] 11.54929
## 
## [[3]]$rect$left
## [1] 11
## 
## [[3]]$rect$top
## [1] 300
## 
## 
## [[3]]$text
## [[3]]$text$x
## [1] 13.05085 13.05085
## 
## [[3]]$text$y
## [1] 296.1502 292.3005
## 
## 
## 
## [[4]]
## [[4]]$rect
## [[4]]$rect$w
## [1] 4.044734
## 
## [[4]]$rect$h
## [1] 36.92899
## 
## [[4]]$rect$left
## [1] 11
## 
## [[4]]$rect$top
## [1] 300
## 
## 
## [[4]]$text
## [[4]]$text$x
## [1] 13.05085 13.05085
## 
## [[4]]$text$y
## [1] 287.6903 275.3807
## 
## 
## 
## [[5]]
## [[5]]$rect
## [[5]]$rect$w
## [1] 4.044734
## 
## [[5]]$rect$h
## [1] 0.8681859
## 
## [[5]]$rect$left
## [1] 11
## 
## [[5]]$rect$top
## [1] 300
## 
## 
## [[5]]$text
## [[5]]$text$x
## [1] 13.05085 13.05085
## 
## [[5]]$text$y
## [1] 299.7106 299.4212
## 
## 
## 
## [[6]]
## [[6]]$rect
## [[6]]$rect$w
## [1] 4.044734
## 
## [[6]]$rect$h
## [1] 0.8347121
## 
## [[6]]$rect$left
## [1] 11
## 
## [[6]]$rect$top
## [1] 300
## 
## 
## [[6]]$text
## [[6]]$text$x
## [1] 13.05085 13.05085
## 
## [[6]]$text$y
## [1] 299.7218 299.4435
## 
## 
## 
## [[7]]
## [[7]]$rect
## [[7]]$rect$w
## [1] 4.044734
## 
## [[7]]$rect$h
## [1] 6.8022
## 
## [[7]]$rect$left
## [1] 11
## 
## [[7]]$rect$top
## [1] 300
## 
## 
## [[7]]$text
## [[7]]$text$x
## [1] 13.05085 13.05085
## 
## [[7]]$text$y
## [1] 297.7326 295.4652
## 
## 
## 
## [[8]]
## [[8]]$rect
## [[8]]$rect$w
## [1] 4.044734
## 
## [[8]]$rect$h
## [1] 3.216695
## 
## [[8]]$rect$left
## [1] 11
## 
## [[8]]$rect$top
## [1] 300
## 
## 
## [[8]]$text
## [[8]]$text$x
## [1] 13.05085 13.05085
## 
## [[8]]$text$y
## [1] 298.9278 297.8555
## 
## 
## 
## [[9]]
## [[9]]$rect
## [[9]]$rect$w
## [1] 4.044734
## 
## [[9]]$rect$h
## [1] 2.940665
## 
## [[9]]$rect$left
## [1] 11
## 
## [[9]]$rect$top
## [1] 300
## 
## 
## [[9]]$text
## [[9]]$text$x
## [1] 13.05085 13.05085
## 
## [[9]]$text$y
## [1] 299.0198 298.0396
## 
## 
## 
## [[10]]
## [[10]]$rect
## [[10]]$rect$w
## [1] 4.044734
## 
## [[10]]$rect$h
## [1] 16.7739
## 
## [[10]]$rect$left
## [1] 11
## 
## [[10]]$rect$top
## [1] 300
## 
## 
## [[10]]$text
## [[10]]$text$x
## [1] 13.05085 13.05085
## 
## [[10]]$text$y
## [1] 294.4087 288.8174
## 
## 
## 
## [[11]]
## [[11]]$rect
## [[11]]$rect$w
## [1] 4.044734
## 
## [[11]]$rect$h
## [1] 2.288299
## 
## [[11]]$rect$left
## [1] 11
## 
## [[11]]$rect$top
## [1] 300
## 
## 
## [[11]]$text
## [[11]]$text$x
## [1] 13.05085 13.05085
## 
## [[11]]$text$y
## [1] 299.2372 298.4745
## 
## 
## 
## [[12]]
## [[12]]$rect
## [[12]]$rect$w
## [1] 4.044734
## 
## [[12]]$rect$h
## [1] 7.829251
## 
## [[12]]$rect$left
## [1] 11
## 
## [[12]]$rect$top
## [1] 300
## 
## 
## [[12]]$text
## [[12]]$text$x
## [1] 13.05085 13.05085
## 
## [[12]]$text$y
## [1] 297.3902 294.7805
## 
## 
## 
## [[13]]
## [[13]]$rect
## [[13]]$rect$w
## [1] 4.044734
## 
## [[13]]$rect$h
## [1] 5.61162
## 
## [[13]]$rect$left
## [1] 11
## 
## [[13]]$rect$top
## [1] 300
## 
## 
## [[13]]$text
## [[13]]$text$x
## [1] 13.05085 13.05085
## 
## [[13]]$text$y
## [1] 298.1295 296.2589
## 
## 
## 
## [[14]]
## [[14]]$rect
## [[14]]$rect$w
## [1] 4.044734
## 
## [[14]]$rect$h
## [1] 2.752561
## 
## [[14]]$rect$left
## [1] 11
## 
## [[14]]$rect$top
## [1] 300
## 
## 
## [[14]]$text
## [[14]]$text$x
## [1] 13.05085 13.05085
## 
## [[14]]$text$y
## [1] 299.0825 298.1650
## 
## 
## 
## [[15]]
## [[15]]$rect
## [[15]]$rect$w
## [1] 4.044734
## 
## [[15]]$rect$h
## [1] 6.043283
## 
## [[15]]$rect$left
## [1] 11
## 
## [[15]]$rect$top
## [1] 300
## 
## 
## [[15]]$text
## [[15]]$text$x
## [1] 13.05085 13.05085
## 
## [[15]]$text$y
## [1] 297.9856 295.9711
## 
## 
## 
## [[16]]
## [[16]]$rect
## [[16]]$rect$w
## [1] 4.044734
## 
## [[16]]$rect$h
## [1] 4.694528
## 
## [[16]]$rect$left
## [1] 11
## 
## [[16]]$rect$top
## [1] 300
## 
## 
## [[16]]$text
## [[16]]$text$x
## [1] 13.05085 13.05085
## 
## [[16]]$text$y
## [1] 298.4352 296.8703
## 
## 
## 
## [[17]]
## [[17]]$rect
## [[17]]$rect$w
## [1] 4.044734
## 
## [[17]]$rect$h
## [1] 1.295886
## 
## [[17]]$rect$left
## [1] 11
## 
## [[17]]$rect$top
## [1] 300
## 
## 
## [[17]]$text
## [[17]]$text$x
## [1] 13.05085 13.05085
## 
## [[17]]$text$y
## [1] 299.5680 299.1361
## 
## 
## 
## [[18]]
## [[18]]$rect
## [[18]]$rect$w
## [1] 4.044734
## 
## [[18]]$rect$h
## [1] 8.090981
## 
## [[18]]$rect$left
## [1] 11
## 
## [[18]]$rect$top
## [1] 300
## 
## 
## [[18]]$text
## [[18]]$text$x
## [1] 13.05085 13.05085
## 
## [[18]]$text$y
## [1] 297.303 294.606
dge2 <- dge[which(rownames(dge) %in% hdacs),]

head(dge2) |>
  kbl(caption = "HDAC genes)") |>
  kable_paper("hover", full_width = F)
HDAC genes)
baseMean log2FoldChange lfcSE stat pvalue padj pos1 pre1 pos2 pre2 pos3 pre3 pos4 pre4 pos5 pre5 pos6 pre6
ENSG00000096717.12 SIRT1 216.46509 0.5519410 0.1484932 3.716944 0.0002016 0.0069943 10.183383 9.976075 10.045656 9.861508 10.112940 10.061779 10.000132 9.853118 10.055607 9.853815 10.108297 9.935662
ENSG00000100429.18 HDAC10 118.93227 -0.5738910 0.2395790 -2.395415 0.0166016 0.1489558 9.770579 9.778175 9.709498 9.772188 9.604090 9.981868 9.940724 10.038299 9.613596 9.829335 9.584435 9.599082
ENSG00000108840.17 HDAC5 929.64390 0.2006982 0.0902478 2.223856 0.0261582 0.1946451 11.044684 11.022195 10.805302 10.710976 10.908164 10.764971 11.238850 10.953372 10.843068 10.868000 10.672234 10.544095
ENSG00000068024.18 HDAC4 796.56869 0.1522836 0.0894729 1.702008 0.0887539 0.3901151 11.019188 11.019509 10.707973 10.717471 10.855695 10.745277 10.666216 10.498403 10.712459 10.596771 10.777771 10.698445
ENSG00000048052.25 HDAC9 380.48088 0.1843314 0.1220415 1.510399 0.1309417 0.4702410 10.168324 10.088455 10.240364 10.305040 10.417715 10.323023 10.477744 10.462694 10.219912 10.086389 10.346098 10.192823
ENSG00000077463.15 SIRT6 48.24593 -0.4012873 0.3278881 -1.223854 0.2210072 0.5906585 9.628205 9.580537 9.499909 9.480289 9.398982 9.510206 9.413109 9.617271 9.634264 9.636789 9.512664 9.587088 
plot(log10(dge$baseMean),dge$log2FoldChange,cex=0.5,pch=19,col="darkgray",
  xlab="log10 basemean",ylab="log2 fold change",main="RNA expression")
points(log10(dge2$baseMean),dge2$log2FoldChange,cex=0.5,pch=19,col="black")
mtext("HDACs")

Session

sessionInfo()
## R version 4.4.3 (2025-02-28)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] network_1.19.0              dplyr_1.1.4                
##  [3] kableExtra_1.4.0            beeswarm_0.4.0             
##  [5] topconfects_1.22.0          limma_3.60.6               
##  [7] eulerr_7.0.2                mitch_1.19.3               
##  [9] gplots_3.2.0                DESeq2_1.44.0              
## [11] SummarizedExperiment_1.34.0 Biobase_2.64.0             
## [13] MatrixGenerics_1.16.0       matrixStats_1.5.0          
## [15] GenomicRanges_1.56.2        GenomeInfoDb_1.40.1        
## [17] IRanges_2.38.1              S4Vectors_0.42.1           
## [19] BiocGenerics_0.50.0         reshape2_1.4.4             
## [21] RhpcBLASctl_0.23-42        
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1        viridisLite_0.4.2       bitops_1.0-9           
##  [4] fastmap_1.2.0           GGally_2.2.1            promises_1.3.2         
##  [7] digest_0.6.37           mime_0.13               lifecycle_1.0.4        
## [10] statmod_1.5.0           magrittr_2.0.3          compiler_4.4.3         
## [13] rlang_1.1.5             sass_0.4.9              tools_4.4.3            
## [16] yaml_2.3.10             knitr_1.50              S4Arrays_1.4.1         
## [19] htmlwidgets_1.6.4       DelayedArray_0.30.1     xml2_1.3.8             
## [22] plyr_1.8.9              RColorBrewer_1.1-3      abind_1.4-8            
## [25] BiocParallel_1.38.0     KernSmooth_2.23-26      purrr_1.0.4            
## [28] grid_4.4.3              caTools_1.18.3          xtable_1.8-4           
## [31] colorspace_2.1-1        ggplot2_3.5.1           MASS_7.3-65            
## [34] scales_1.3.0            gtools_3.9.5            cli_3.6.4              
## [37] rmarkdown_2.29          crayon_1.5.3            generics_0.1.3         
## [40] rstudioapi_0.17.1       httr_1.4.7              cachem_1.1.0           
## [43] stringr_1.5.1           zlibbioc_1.50.0         assertthat_0.2.1       
## [46] parallel_4.4.3          XVector_0.44.0          vctrs_0.6.5            
## [49] Matrix_1.7-3            jsonlite_1.9.1          echarts4r_0.4.5        
## [52] systemfonts_1.2.1       locfit_1.5-9.12         jquerylib_0.1.4        
## [55] tidyr_1.3.1             glue_1.8.0              statnet.common_4.11.0  
## [58] ggstats_0.9.0           codetools_0.2-20        stringi_1.8.4          
## [61] gtable_0.3.6            later_1.4.1             UCSC.utils_1.0.0       
## [64] munsell_0.5.1           tibble_3.2.1            pillar_1.10.1          
## [67] htmltools_0.5.8.1       GenomeInfoDbData_1.2.12 R6_2.6.1               
## [70] evaluate_1.0.3          shiny_1.10.0            lattice_0.22-6         
## [73] httpuv_1.6.15           bslib_0.9.0             Rcpp_1.0.14            
## [76] coda_0.19-4.1           svglite_2.1.3           gridExtra_2.3          
## [79] SparseArray_1.4.8       xfun_0.51               pkgconfig_2.0.3