Introduction

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

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)

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)

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

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

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)

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

Session

sessionInfo()
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 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] kableExtra_1.4.0            beeswarm_0.4.0             
##  [3] topconfects_1.20.0          limma_3.60.0               
##  [5] eulerr_7.0.2                mitch_1.16.0               
##  [7] gplots_3.1.3.1              DESeq2_1.44.0              
##  [9] SummarizedExperiment_1.34.0 Biobase_2.64.0             
## [11] MatrixGenerics_1.16.0       matrixStats_1.3.0          
## [13] GenomicRanges_1.56.0        GenomeInfoDb_1.40.0        
## [15] IRanges_2.38.0              S4Vectors_0.42.0           
## [17] BiocGenerics_0.50.0         reshape2_1.4.4             
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1        viridisLite_0.4.2       dplyr_1.1.4            
##  [4] bitops_1.0-7            fastmap_1.2.0           GGally_2.2.1           
##  [7] promises_1.3.0          digest_0.6.36           mime_0.12              
## [10] lifecycle_1.0.4         statmod_1.5.0           magrittr_2.0.3         
## [13] compiler_4.4.1          rlang_1.1.4             sass_0.4.9             
## [16] tools_4.4.1             utf8_1.2.4              yaml_2.3.9             
## [19] knitr_1.48              S4Arrays_1.4.0          htmlwidgets_1.6.4      
## [22] DelayedArray_0.30.1     xml2_1.3.6              plyr_1.8.9             
## [25] RColorBrewer_1.1-3      abind_1.4-5             BiocParallel_1.38.0    
## [28] KernSmooth_2.23-24      purrr_1.0.2             grid_4.4.1             
## [31] fansi_1.0.6             caTools_1.18.2          xtable_1.8-4           
## [34] colorspace_2.1-0        ggplot2_3.5.1           MASS_7.3-61            
## [37] scales_1.3.0            gtools_3.9.5            cli_3.6.3              
## [40] rmarkdown_2.27          crayon_1.5.3            generics_0.1.3         
## [43] rstudioapi_0.16.0       httr_1.4.7              cachem_1.1.0           
## [46] stringr_1.5.1           zlibbioc_1.50.0         assertthat_0.2.1       
## [49] parallel_4.4.1          XVector_0.44.0          vctrs_0.6.5            
## [52] Matrix_1.7-0            jsonlite_1.8.8          echarts4r_0.4.5        
## [55] systemfonts_1.1.0       locfit_1.5-9.10         jquerylib_0.1.4        
## [58] tidyr_1.3.1             glue_1.7.0              ggstats_0.6.0          
## [61] codetools_0.2-20        stringi_1.8.4           gtable_0.3.5           
## [64] later_1.3.2             UCSC.utils_1.0.0        munsell_0.5.1          
## [67] tibble_3.2.1            pillar_1.9.0            htmltools_0.5.8.1      
## [70] GenomeInfoDbData_1.2.12 R6_2.5.1                evaluate_0.24.0        
## [73] shiny_1.8.1.1           lattice_0.22-6          highr_0.11             
## [76] httpuv_1.6.15           bslib_0.7.0             Rcpp_1.0.12            
## [79] svglite_2.1.3           gridExtra_2.3           SparseArray_1.4.3      
## [82] xfun_0.45               pkgconfig_2.0.3