Exercise RNA-seq
Mark Ziemann
2025-02-21
Introduction
suppressPackageStartupMessages({
library("RhpcBLASctl")
library("reshape2")
library("DESeq2")
library("gplots")
library("mitch")
library("eulerr")
library("limma")
library("topconfects")
library("beeswarm")
library("kableExtra")
})
RhpcBLASctl::blas_set_num_threads(1)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 12colSums(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 20436123colnames(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()## X11cairo
## 2DESeq2
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 12xx <- xx[rowMeans(xx)>=10,]
dim(xx)## [1] 16971 12dds <- DESeqDataSetFromMatrix(countData = xx , colData = ss, design = ~ pos )## converting counts to integer moderes <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- 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)| 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] 129nrow(subset(dge_unpaired,padj<0.05 & log2FoldChange < 0))## [1] 28dds <- 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 factorsres <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingz <- 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)| 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] 472nrow(subset(dge,padj<0.05 & log2FoldChange < 0))## [1] 165Smear 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()## X11cairo
## 2Mitch
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.997head(m)## x
## 5_8S_rRNA 2.8723535
## 7SK 0.1300210
## A1BG 0.2484592
## A2M 3.4477904
## A2M-AS1 -0.5669178
## A4GALT 0.1378113mres <- 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)| 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)| 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] 2373up <- 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()## X11cairo
## 2par(mar=c(5.1, 4.1, 4.1, 2.1))
if (! file.exists("rna_mitch.html") ) {
mitch_report(mres,"rna_mitch.html")
}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.2494SLC2A4 <- 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()## X11cairo
## 2tpm <- apply(x[,1:12], 2 , function(x) { x / sum(x) } ) * 1000000
tpm <- tpm[rowMeans(tpm) > 1,]
dim(tpm)## [1] 36942 12SLC2A4 <- 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.16995HATs 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.8035dge1 <- dge[which(rownames(dge) %in% hats),]
head(dge1) |>
kbl(caption = "HAT genes)") |>
kable_paper("hover", full_width = F)| 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.606dge2 <- dge[which(rownames(dge) %in% hdacs),]
head(dge2) |>
kbl(caption = "HDAC genes)") |>
kable_paper("hover", full_width = F)| 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.2 (2024-10-31)
## 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] pkgload_1.4.0 GGally_2.2.1
## [3] ggplot2_3.5.1 gtools_3.9.5
## [5] tibble_3.2.1 dplyr_1.1.4
## [7] echarts4r_0.4.5 kableExtra_1.4.0
## [9] beeswarm_0.4.0 topconfects_1.22.0
## [11] limma_3.60.6 eulerr_7.0.2
## [13] mitch_1.19.3 gplots_3.2.0
## [15] DESeq2_1.44.0 SummarizedExperiment_1.34.0
## [17] Biobase_2.64.0 MatrixGenerics_1.16.0
## [19] matrixStats_1.4.1 GenomicRanges_1.56.2
## [21] GenomeInfoDb_1.40.1 IRanges_2.38.1
## [23] S4Vectors_0.42.1 BiocGenerics_0.50.0
## [25] reshape2_1.4.4 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 promises_1.3.1 digest_0.6.37
## [7] mime_0.12 lifecycle_1.0.4 statmod_1.5.0
## [10] magrittr_2.0.3 compiler_4.4.2 rlang_1.1.4
## [13] sass_0.4.9 tools_4.4.2 utf8_1.2.4
## [16] yaml_2.3.10 knitr_1.49 S4Arrays_1.4.1
## [19] htmlwidgets_1.6.4 DelayedArray_0.30.1 xml2_1.3.6
## [22] plyr_1.8.9 RColorBrewer_1.1-3 abind_1.4-8
## [25] BiocParallel_1.38.0 KernSmooth_2.23-24 withr_3.0.2
## [28] purrr_1.0.2 grid_4.4.2 fansi_1.0.6
## [31] caTools_1.18.3 xtable_1.8-4 colorspace_2.1-1
## [34] MASS_7.3-64 scales_1.3.0 cli_3.6.3
## [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.2 XVector_0.44.0 vctrs_0.6.5
## [49] Matrix_1.7-1 jsonlite_1.8.9 systemfonts_1.1.0
## [52] locfit_1.5-9.10 jquerylib_0.1.4 tidyr_1.3.1
## [55] glue_1.8.0 ggstats_0.7.0 codetools_0.2-20
## [58] stringi_1.8.4 gtable_0.3.6 later_1.4.0
## [61] UCSC.utils_1.0.0 munsell_0.5.1 pillar_1.9.0
## [64] htmltools_0.5.8.1 GenomeInfoDbData_1.2.12 R6_2.5.1
## [67] evaluate_1.0.1 shiny_1.9.1 lattice_0.22-6
## [70] httpuv_1.6.15 bslib_0.8.0 Rcpp_1.0.13-1
## [73] svglite_2.1.3 gridExtra_2.3 SparseArray_1.4.8
## [76] xfun_0.49 pkgconfig_2.0.3