Exercise ChIP-seq
Mark Ziemann
2024-12-02
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 and analyse read positions
Firstly for H3K9ac.
myfiles <- list.files("raw_data/h3k9ac/",pattern="*pos2",full.names=TRUE)
res <- sapply(myfiles,function(f) {
x <- read.table(f,header=FALSE)
x <- x[grep("chr",x$V1),]
x <- subset(x, V8 < 5000 & V8 > -5000 )
xh <- hist(x$V8,breaks=200)
xh$counts
} )
inputs <- res[,grep("SRR",colnames(res))]
res <- res[,grep("SRR",colnames(res),invert=TRUE)]
# normalise
y <- apply(res , 2, function(x) { x/sum(x)*100 } )
# plot all
plot(-100:99*50,res[,1]/sum(res[,1])*100,type="p",cex=0,ylim=c(0,3),
xlab="Distance to TSS", ylab="read density")
lapply(1:ncol(res) , function(i) {
y <- res[,i]/sum(res[,i])*100
points(-100:99*50,y,type="l",col=i%%2+1)
} )## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
##
## [[7]]
## NULL
##
## [[8]]
## NULL
##
## [[9]]
## NULL
##
## [[10]]
## NULL
##
## [[11]]
## NULL
##
## [[12]]
## NULLlapply(1:ncol(inputs) , function(i) {
y <- inputs[,i]/sum(inputs[,i])*100
points(-100:99*50,y,type="l",col="blue")
} )
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL# plot each
par(mfrow=c(3,2))
plot(-100:99*50,y[,2],type="l",cex=1,ylim=c(0,3), main="P1",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*50,y[,1],type="l",col=2)
plot(-100:99*50,y[,4],type="l",cex=1,ylim=c(0,3), main="P2",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*50,y[,3],type="l",col=2)
plot(-100:99*50,y[,6],type="l",cex=1,ylim=c(0,3), main="P3",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*50,y[,5],type="l",col=2)
plot(-100:99*50,y[,8],type="l",cex=1,ylim=c(0,3), main="P4",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*50,y[,7],type="l",col=2)
plot(-100:99*50,y[,10],type="l",cex=1,ylim=c(0,3), main="P5",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*50,y[,9],type="l",col=2)
plot(-100:99*50,y[,12],type="l",cex=1,ylim=c(0,3), main="P6",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*50,y[,11],type="l",col=2)
par(mfrow=c(1,1))
k9pre <- rowMeans(y[,seq(2,12,2)])
k9pos <- rowMeans(y[,seq(1,12,2)])Now for H3K36ac.
myfiles <- list.files("raw_data/h3k36ac/",pattern="*pos2",full.names=TRUE)
res <- sapply(myfiles,function(f) {
x <- read.table(f,header=FALSE)
x <- x[grep("chr",x$V1),]
x <- subset(x, V8 < 5000 & V8 > -5000 )
xh <- hist(x$V8,breaks=200)
xh$counts
} )
# normalise
y <- apply(res , 2, function(x) { x/sum(x)*100 } )
# plot all
plot(-100:99*50,res[,1]/sum(res[,1])*100,type="p",cex=0,ylim=c(0.3,1.7),
xlab="Distance to TSS", ylab="read density")
lapply(1:ncol(res) , function(i) {
y <- res[,i]/sum(res[,i])*100
points(-100:99*50,y,type="l",col=i%%2+1)
} )
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULL
##
## [[5]]
## NULL
##
## [[6]]
## NULL
##
## [[7]]
## NULL
##
## [[8]]
## NULL
##
## [[9]]
## NULL
##
## [[10]]
## NULL
##
## [[11]]
## NULL
##
## [[12]]
## NULL# plot each
par(mfrow=c(3,2))
plot(-100:99*10,y[,2],type="l",cex=1,ylim=c(0.3,1.7), main="P1",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*10,y[,1],type="l",col=2)
plot(-100:99*10,y[,4],type="l",cex=1,ylim=c(0.3,1.7), main="P2",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*10,y[,3],type="l",col=2)
plot(-100:99*10,y[,6],type="l",cex=1,ylim=c(0.3,1.7), main="P3",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*10,y[,5],type="l",col=2)
plot(-100:99*10,y[,8],type="l",cex=1,ylim=c(0.3,1.7), main="P4",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*10,y[,7],type="l",col=2)
plot(-100:99*10,y[,10],type="l",cex=1,ylim=c(0.3,1.7), main="P5",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*10,y[,9],type="l",col=2)
plot(-100:99*10,y[,12],type="l",cex=1,ylim=c(0.3,1.7), main="P6",
xlab="Distance to TSS", ylab="read density",col=4)
points(-100:99*10,y[,11],type="l",col=2)
k36pre <- rowMeans(y[,seq(2,12,2)])
k36pos <- rowMeans(y[,seq(1,12,2)])
par(mfrow=c(1,1))Now plot mean pre for both k9 and k36.
par(mfrow=c(2,1))
plot(-100:99*50,k9pre,type="l",cex=1,ylim=c(0.3,2.2),bty="n",lwd=2,
xlab="Distance to TSS", ylab="read density",col=4)
mtext("H3K9ac")
points(-100:99*50,k9pos,type="l",col=2,lwd=2)
legend("topright", legend=c("Post", "Pre"),
col=c(2,4), lty=1, cex=0.9, lwd=2)
plot(-100:99*50,k36pre,type="l",cex=1,ylim=c(0.3,1.5),bty="n",lwd=2,
xlab="Distance to TSS", ylab="read density",col=4)
mtext("H3K36ac")
points(-100:99*50,k36pos,type="l",col=2,lwd=2)
legend("topright", legend=c("Post", "Pre"),
col=c(2,4), lty=1, cex=0.9, lwd=2)
par(mfrow=c(1,1))
pdf("chip_tss_mds.pdf",width=5,height=5)
par(mfrow=c(2,1))
plot(-100:99*50,k9pre,type="l",cex=1,ylim=c(0.3,2.2),bty="n",lwd=2,
xlab="Distance to TSS", ylab="read density",col=4)
mtext("H3K9ac")
points(-100:99*50,k9pos,type="l",col=2,lwd=2)
legend("topright", legend=c("Post", "Pre"),
col=c(2,4), lty=1, cex=0.9, lwd=2)
plot(-100:99*50,k36pre,type="l",cex=1,ylim=c(0.3,1.5),bty="n",lwd=2,
xlab="Distance to TSS", ylab="read density",col=4)
mtext("H3K36ac")
points(-100:99*50,k36pos,type="l",col=2,lwd=2)
legend("topright", legend=c("Post", "Pre"),
col=c(2,4), lty=1, cex=0.9, lwd=2)
par(mfrow=c(1,1))
dev.off()## png
## 2Import read counts
Rename samples. Don’t use initials.
Get the proportion of reads that are within 2kbp of TSSs.
k9s <- read.table("raw_data/h3k9ac/h3k9_counts.tsv.summary",header=TRUE,row.names=1)
inputs <- k9s[,grep("SRR",colnames(k9s))]
k9s <- k9s[,grep("SRR",colnames(k9s),invert=TRUE)]
colnames(k9s) <- c("pos1","pre1","pos2","pre2","pos3","pre3","pos4","pre4","pos5","pre5","pos6","pre6")
colnames(inputs) <- c("input1","input2","input3","input4")
inputs## input1 input2 input3 input4
## Assigned 3888522 4025908 2245486 2159094
## Unassigned_Unmapped 253745 222771 452959 2863289
## Unassigned_Read_Type 0 0 0 0
## Unassigned_Singleton 0 0 0 0
## Unassigned_MappingQuality 2879803 2551110 1894829 1564142
## Unassigned_Chimera 0 0 0 0
## Unassigned_FragmentLength 0 0 0 0
## Unassigned_Duplicate 0 0 0 0
## Unassigned_MultiMapping 0 0 0 0
## Unassigned_Secondary 0 0 0 0
## Unassigned_NonSplit 0 0 0 0
## Unassigned_NoFeatures 21620001 21494794 12686447 12199144
## Unassigned_Overlapping_Length 0 0 0 0
## Unassigned_Ambiguity 0 0 0 0barplot(as.matrix(inputs))
inputsn <- apply(inputs,2,function(x) {x/sum(x)} )
inputsn## input1 input2 input3 input4
## Assigned 0.135762599 0.142285469 0.12994920 0.11493304
## Unassigned_Unmapped 0.008859171 0.007873274 0.02621333 0.15241879
## Unassigned_Read_Type 0.000000000 0.000000000 0.00000000 0.00000000
## Unassigned_Singleton 0.000000000 0.000000000 0.00000000 0.00000000
## Unassigned_MappingQuality 0.100544510 0.090162488 0.10965623 0.08326251
## Unassigned_Chimera 0.000000000 0.000000000 0.00000000 0.00000000
## Unassigned_FragmentLength 0.000000000 0.000000000 0.00000000 0.00000000
## Unassigned_Duplicate 0.000000000 0.000000000 0.00000000 0.00000000
## Unassigned_MultiMapping 0.000000000 0.000000000 0.00000000 0.00000000
## Unassigned_Secondary 0.000000000 0.000000000 0.00000000 0.00000000
## Unassigned_NonSplit 0.000000000 0.000000000 0.00000000 0.00000000
## Unassigned_NoFeatures 0.754833720 0.759678770 0.73418124 0.64938566
## Unassigned_Overlapping_Length 0.000000000 0.000000000 0.00000000 0.00000000
## Unassigned_Ambiguity 0.000000000 0.000000000 0.00000000 0.00000000barplot(as.matrix(inputsn))
k9s## pos1 pre1 pos2 pre2 pos3
## Assigned 7756749 6430759 13082082 4953269 7726068
## Unassigned_Unmapped 385052 557590 751361 585880 412283
## Unassigned_Read_Type 0 0 0 0 0
## Unassigned_Singleton 0 0 0 0 0
## Unassigned_MappingQuality 5725774 7754993 11636911 6738231 6988437
## Unassigned_Chimera 0 0 0 0 0
## Unassigned_FragmentLength 0 0 0 0 0
## Unassigned_Duplicate 0 0 0 0 0
## Unassigned_MultiMapping 0 0 0 0 0
## Unassigned_Secondary 0 0 0 0 0
## Unassigned_NonSplit 0 0 0 0 0
## Unassigned_NoFeatures 29812445 38616775 58178138 30628975 34489472
## Unassigned_Overlapping_Length 0 0 0 0 0
## Unassigned_Ambiguity 0 0 0 0 0
## pre3 pos4 pre4 pos5 pre5
## Assigned 5400705 7540527 11724157 6200147 6033734
## Unassigned_Unmapped 572918 390416 977005 496997 1138925
## Unassigned_Read_Type 0 0 0 0 0
## Unassigned_Singleton 0 0 0 0 0
## Unassigned_MappingQuality 6378033 6532301 13351773 7041413 6293562
## Unassigned_Chimera 0 0 0 0 0
## Unassigned_FragmentLength 0 0 0 0 0
## Unassigned_Duplicate 0 0 0 0 0
## Unassigned_MultiMapping 0 0 0 0 0
## Unassigned_Secondary 0 0 0 0 0
## Unassigned_NonSplit 0 0 0 0 0
## Unassigned_NoFeatures 31630276 33917230 61155197 35344376 30857841
## Unassigned_Overlapping_Length 0 0 0 0 0
## Unassigned_Ambiguity 0 0 0 0 0
## pos6 pre6
## Assigned 8137338 15668690
## Unassigned_Unmapped 747681 1918477
## Unassigned_Read_Type 0 0
## Unassigned_Singleton 0 0
## Unassigned_MappingQuality 7006344 16463353
## Unassigned_Chimera 0 0
## Unassigned_FragmentLength 0 0
## Unassigned_Duplicate 0 0
## Unassigned_MultiMapping 0 0
## Unassigned_Secondary 0 0
## Unassigned_NonSplit 0 0
## Unassigned_NoFeatures 34904688 73176720
## Unassigned_Overlapping_Length 0 0
## Unassigned_Ambiguity 0 0barplot(as.matrix(k9s))
k9sn <- apply(k9s,2,function(x) {x/sum(x)} )
k9sn## pos1 pre1 pos2 pre2
## Assigned 0.177581169 0.12051621 0.156393519 0.11544371
## Unassigned_Unmapped 0.008815289 0.01044956 0.008982362 0.01365485
## Unassigned_Read_Type 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_Singleton 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_MappingQuality 0.131084510 0.14533313 0.139116806 0.15704506
## Unassigned_Chimera 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_FragmentLength 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_Duplicate 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_MultiMapping 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_Secondary 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_NonSplit 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_NoFeatures 0.682519033 0.72370109 0.695507314 0.71385637
## Unassigned_Overlapping_Length 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_Ambiguity 0.000000000 0.00000000 0.000000000 0.00000000
## pos3 pre3 pos4 pre4
## Assigned 0.155716453 0.12279372 0.155858891 0.13443880
## Unassigned_Unmapped 0.008309433 0.01302621 0.008069702 0.01120314
## Unassigned_Read_Type 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_Singleton 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_MappingQuality 0.140849734 0.14501484 0.135019368 0.15310238
## Unassigned_Chimera 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_FragmentLength 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_Duplicate 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_MultiMapping 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_Secondary 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_NonSplit 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_NoFeatures 0.695124381 0.71916522 0.701052040 0.70125567
## Unassigned_Overlapping_Length 0.000000000 0.00000000 0.000000000 0.00000000
## Unassigned_Ambiguity 0.000000000 0.00000000 0.000000000 0.00000000
## pos5 pre5 pos6 pre6
## Assigned 0.12631981 0.13612773 0.16019627 0.14612602
## Unassigned_Unmapped 0.01012566 0.02569541 0.01471927 0.01789169
## Unassigned_Read_Type 0.00000000 0.00000000 0.00000000 0.00000000
## Unassigned_Singleton 0.00000000 0.00000000 0.00000000 0.00000000
## Unassigned_MappingQuality 0.14345950 0.14198974 0.13793088 0.15353704
## Unassigned_Chimera 0.00000000 0.00000000 0.00000000 0.00000000
## Unassigned_FragmentLength 0.00000000 0.00000000 0.00000000 0.00000000
## Unassigned_Duplicate 0.00000000 0.00000000 0.00000000 0.00000000
## Unassigned_MultiMapping 0.00000000 0.00000000 0.00000000 0.00000000
## Unassigned_Secondary 0.00000000 0.00000000 0.00000000 0.00000000
## Unassigned_NonSplit 0.00000000 0.00000000 0.00000000 0.00000000
## Unassigned_NoFeatures 0.72009503 0.69618712 0.68715357 0.68244524
## Unassigned_Overlapping_Length 0.00000000 0.00000000 0.00000000 0.00000000
## Unassigned_Ambiguity 0.00000000 0.00000000 0.00000000 0.00000000barplot(as.matrix(k9sn))
k36s <- read.table("raw_data/h3k36ac/h3k36_counts.tsv.summary",header=TRUE,row.names=1)
colnames(k36s) <- c("pos1","pre1","pos2","pre2","pos3","pre3","pos4","pre4","pos5","pre5","pos6","pre6")
k36s## pos1 pre1 pos2 pre2 pos3
## Assigned 6195544 6294768 5618086 5303142 5733831
## Unassigned_Unmapped 1473465 1111917 700103 841680 1774859
## Unassigned_Read_Type 0 0 0 0 0
## Unassigned_Singleton 0 0 0 0 0
## Unassigned_MappingQuality 9524597 9590123 9109908 8233487 8441755
## Unassigned_Chimera 0 0 0 0 0
## Unassigned_FragmentLength 0 0 0 0 0
## Unassigned_Duplicate 0 0 0 0 0
## Unassigned_MultiMapping 0 0 0 0 0
## Unassigned_Secondary 0 0 0 0 0
## Unassigned_NonSplit 0 0 0 0 0
## Unassigned_NoFeatures 42129671 40474063 40097125 35343506 36537911
## Unassigned_Overlapping_Length 0 0 0 0 0
## Unassigned_Ambiguity 0 0 0 0 0
## pre3 pos4 pre4 pos5 pre5
## Assigned 6664455 6926829 8821544 6137079 6696549
## Unassigned_Unmapped 1556532 2016880 935526 918295 744888
## Unassigned_Read_Type 0 0 0 0 0
## Unassigned_Singleton 0 0 0 0 0
## Unassigned_MappingQuality 10367327 10774766 12242494 9369415 8997966
## Unassigned_Chimera 0 0 0 0 0
## Unassigned_FragmentLength 0 0 0 0 0
## Unassigned_Duplicate 0 0 0 0 0
## Unassigned_MultiMapping 0 0 0 0 0
## Unassigned_Secondary 0 0 0 0 0
## Unassigned_NonSplit 0 0 0 0 0
## Unassigned_NoFeatures 43655931 48270207 46488183 41198373 33778390
## Unassigned_Overlapping_Length 0 0 0 0 0
## Unassigned_Ambiguity 0 0 0 0 0
## pos6 pre6
## Assigned 5823145 6565429
## Unassigned_Unmapped 912349 1071548
## Unassigned_Read_Type 0 0
## Unassigned_Singleton 0 0
## Unassigned_MappingQuality 8906465 8984204
## Unassigned_Chimera 0 0
## Unassigned_FragmentLength 0 0
## Unassigned_Duplicate 0 0
## Unassigned_MultiMapping 0 0
## Unassigned_Secondary 0 0
## Unassigned_NonSplit 0 0
## Unassigned_NoFeatures 38244027 34112607
## Unassigned_Overlapping_Length 0 0
## Unassigned_Ambiguity 0 0barplot(as.matrix(k36s))
k36sn <- apply(k36s,2,function(x) {x/sum(x)} )
barplot(as.matrix(k36sn))
all <- cbind(inputs,k9s,k36s)
all## input1 input2 input3 input4 pos1
## Assigned 3888522 4025908 2245486 2159094 7756749
## Unassigned_Unmapped 253745 222771 452959 2863289 385052
## Unassigned_Read_Type 0 0 0 0 0
## Unassigned_Singleton 0 0 0 0 0
## Unassigned_MappingQuality 2879803 2551110 1894829 1564142 5725774
## Unassigned_Chimera 0 0 0 0 0
## Unassigned_FragmentLength 0 0 0 0 0
## Unassigned_Duplicate 0 0 0 0 0
## Unassigned_MultiMapping 0 0 0 0 0
## Unassigned_Secondary 0 0 0 0 0
## Unassigned_NonSplit 0 0 0 0 0
## Unassigned_NoFeatures 21620001 21494794 12686447 12199144 29812445
## Unassigned_Overlapping_Length 0 0 0 0 0
## Unassigned_Ambiguity 0 0 0 0 0
## pre1 pos2 pre2 pos3 pre3
## Assigned 6430759 13082082 4953269 7726068 5400705
## Unassigned_Unmapped 557590 751361 585880 412283 572918
## Unassigned_Read_Type 0 0 0 0 0
## Unassigned_Singleton 0 0 0 0 0
## Unassigned_MappingQuality 7754993 11636911 6738231 6988437 6378033
## Unassigned_Chimera 0 0 0 0 0
## Unassigned_FragmentLength 0 0 0 0 0
## Unassigned_Duplicate 0 0 0 0 0
## Unassigned_MultiMapping 0 0 0 0 0
## Unassigned_Secondary 0 0 0 0 0
## Unassigned_NonSplit 0 0 0 0 0
## Unassigned_NoFeatures 38616775 58178138 30628975 34489472 31630276
## Unassigned_Overlapping_Length 0 0 0 0 0
## Unassigned_Ambiguity 0 0 0 0 0
## pos4 pre4 pos5 pre5 pos6
## Assigned 7540527 11724157 6200147 6033734 8137338
## Unassigned_Unmapped 390416 977005 496997 1138925 747681
## Unassigned_Read_Type 0 0 0 0 0
## Unassigned_Singleton 0 0 0 0 0
## Unassigned_MappingQuality 6532301 13351773 7041413 6293562 7006344
## Unassigned_Chimera 0 0 0 0 0
## Unassigned_FragmentLength 0 0 0 0 0
## Unassigned_Duplicate 0 0 0 0 0
## Unassigned_MultiMapping 0 0 0 0 0
## Unassigned_Secondary 0 0 0 0 0
## Unassigned_NonSplit 0 0 0 0 0
## Unassigned_NoFeatures 33917230 61155197 35344376 30857841 34904688
## Unassigned_Overlapping_Length 0 0 0 0 0
## Unassigned_Ambiguity 0 0 0 0 0
## pre6 pos1 pre1 pos2 pre2
## Assigned 15668690 6195544 6294768 5618086 5303142
## Unassigned_Unmapped 1918477 1473465 1111917 700103 841680
## Unassigned_Read_Type 0 0 0 0 0
## Unassigned_Singleton 0 0 0 0 0
## Unassigned_MappingQuality 16463353 9524597 9590123 9109908 8233487
## Unassigned_Chimera 0 0 0 0 0
## Unassigned_FragmentLength 0 0 0 0 0
## Unassigned_Duplicate 0 0 0 0 0
## Unassigned_MultiMapping 0 0 0 0 0
## Unassigned_Secondary 0 0 0 0 0
## Unassigned_NonSplit 0 0 0 0 0
## Unassigned_NoFeatures 73176720 42129671 40474063 40097125 35343506
## Unassigned_Overlapping_Length 0 0 0 0 0
## Unassigned_Ambiguity 0 0 0 0 0
## pos3 pre3 pos4 pre4 pos5
## Assigned 5733831 6664455 6926829 8821544 6137079
## Unassigned_Unmapped 1774859 1556532 2016880 935526 918295
## Unassigned_Read_Type 0 0 0 0 0
## Unassigned_Singleton 0 0 0 0 0
## Unassigned_MappingQuality 8441755 10367327 10774766 12242494 9369415
## Unassigned_Chimera 0 0 0 0 0
## Unassigned_FragmentLength 0 0 0 0 0
## Unassigned_Duplicate 0 0 0 0 0
## Unassigned_MultiMapping 0 0 0 0 0
## Unassigned_Secondary 0 0 0 0 0
## Unassigned_NonSplit 0 0 0 0 0
## Unassigned_NoFeatures 36537911 43655931 48270207 46488183 41198373
## Unassigned_Overlapping_Length 0 0 0 0 0
## Unassigned_Ambiguity 0 0 0 0 0
## pre5 pos6 pre6
## Assigned 6696549 5823145 6565429
## Unassigned_Unmapped 744888 912349 1071548
## Unassigned_Read_Type 0 0 0
## Unassigned_Singleton 0 0 0
## Unassigned_MappingQuality 8997966 8906465 8984204
## Unassigned_Chimera 0 0 0
## Unassigned_FragmentLength 0 0 0
## Unassigned_Duplicate 0 0 0
## Unassigned_MultiMapping 0 0 0
## Unassigned_Secondary 0 0 0
## Unassigned_NonSplit 0 0 0
## Unassigned_NoFeatures 33778390 38244027 34112607
## Unassigned_Overlapping_Length 0 0 0
## Unassigned_Ambiguity 0 0 0allsn <- apply(all,2,function(x) {x/sum(x)} )
barplot(as.matrix(allsn),ylab="proportion of reads",las=2)
k9tss <- k9sn[1,] / (k9sn[1,] + k9sn[12,] )
inputtss <- unlist(as.vector(inputs[1,] / ( inputs[1,] + inputs[12,] )))
k36tss <- k36sn[1,] / (k36sn[1,] + k36sn[12,] )
k9tss_pre <- k9tss[grep("pre",names(k9tss))]
k9tss_post <- k9tss[grep("pos",names(k9tss))]
k36tss_pre <- k36tss[grep("pre",names(k36tss))]
k36tss_post <- k36tss[grep("pos",names(k36tss))]
xl <- list("Input"=inputtss,"K9ac pre"=k9tss_pre, "K9ac post"=k9tss_post,
"K36ac pre"=k36tss_pre, "K36ac post"=k36tss_post)
boxplot(xl,cex=0,ylab="Proportion",main="TSS read fraction",col="white")
beeswarm(xl,add=TRUE,pch=19,cex=1.5,col="darkgray")
pdf("chip_boxplot_mds.pdf",width=5,height=5)
boxplot(xl,cex=0,ylab="Proportion",main="TSS read fraction",col="white")
beeswarm(xl,add=TRUE,pch=19,cex=1.5,col="darkgray")
dev.off()## png
## 2Inport data.
gt <- read.table("gt.tsv")
gt$ens <- sapply(strsplit(gt[,1]," "),"[[",1)
dim(gt)## [1] 63086 3k9 <- read.table("raw_data/h3k9ac/h3k9_counts.tsv",header=TRUE)
dim(k9)## [1] 95874 22k9a <- merge(k9,gt,by.x="Geneid",by.y="ens",all.x=TRUE)
rownames(k9a) <- paste(k9a$Geneid,k9a$genesymbol,k9a$Chr,k9a$Start,k9a$End)
k9a <- k9a[,7:22]
colnames(k9a) <- c("pos1","pre1","pos2","pre2","pos3","pre3","pos4","pre4","pos5","pre5",
"input1","input2","input3","input4","pos6","pre6")
input <- k9a[,grep("input",colnames(k9a))]
k9a <- k9a[,grep("input",colnames(k9a),invert=TRUE)]
head(k9a) |>
kbl(caption = "H3K9ac TSS data") |>
kable_paper("hover", full_width = F)| pos1 | pre1 | pos2 | pre2 | pos3 | pre3 | pos4 | pre4 | pos5 | pre5 | pos6 | pre6 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ENSG00000000003.16 TSPAN6 chrX 100634689 100641991 | 107 | 76 | 207 | 57 | 101 | 70 | 102 | 147 | 73 | 114 | 129 | 311 |
| ENSG00000000005.6 TNMD chrX 100582936 100587019 | 32 | 37 | 71 | 34 | 43 | 40 | 39 | 74 | 52 | 41 | 33 | 102 |
| ENSG00000000005.6 TNMD chrX 100591624 100595624 | 30 | 18 | 69 | 33 | 36 | 31 | 38 | 74 | 36 | 21 | 22 | 11 |
| ENSG00000000419.14 DPM1 chr20 50940869 50951472 | 171 | 141 | 275 | 68 | 142 | 85 | 142 | 181 | 112 | 106 | 142 | 239 |
| ENSG00000000419.14 DPM1 chr20 50955153 50961140 | 389 | 140 | 545 | 123 | 283 | 137 | 342 | 307 | 260 | 223 | 343 | 509 |
| ENSG00000000457.14 SCYL3 chr1 169891896 169896267 | 281 | 158 | 444 | 113 | 293 | 146 | 251 | 349 | 149 | 153 | 330 | 368 |
head(input) |>
kbl(caption = "Input TSS data") |>
kable_paper("hover", full_width = F)| input1 | input2 | input3 | input4 | |
|---|---|---|---|---|
| ENSG00000000003.16 TSPAN6 chrX 100634689 100641991 | 74 | 81 | 22 | 30 |
| ENSG00000000005.6 TNMD chrX 100582936 100587019 | 50 | 40 | 11 | 5 |
| ENSG00000000005.6 TNMD chrX 100591624 100595624 | 25 | 25 | 8 | 0 |
| ENSG00000000419.14 DPM1 chr20 50940869 50951472 | 138 | 121 | 72 | 81 |
| ENSG00000000419.14 DPM1 chr20 50955153 50961140 | 97 | 93 | 31 | 62 |
| ENSG00000000457.14 SCYL3 chr1 169891896 169896267 | 66 | 65 | 40 | 18 |
colSums(k9a)## pos1 pre1 pos2 pre2 pos3 pre3 pos4 pre4
## 9922023 7718857 16458895 5881785 9807888 6507463 9553319 14483957
## pos5 pre5 pos6 pre6
## 7508858 7411365 10334114 19855186k36 <- read.table("raw_data/h3k36ac/h3k36_counts.tsv",header=TRUE)
dim(k36)## [1] 95874 18k36a <- merge(k36,gt,by.x="Geneid",by.y="ens",all.x=TRUE)
rownames(k36a) <- paste(k36a$Geneid,k36a$genesymbol,k36a$Chr,k36a$Start,k36a$End)
k36a <- k36a[,7:18]
colnames(k36a) <- c("pos1","pre1","pos2","pre2","pos3","pre3","pos4","pre4","pos5","pre5","pos6","pre6")
head(k36a) |>
kbl(caption = "H3K36ac TSS data") |>
kable_paper("hover", full_width = F)| pos1 | pre1 | pos2 | pre2 | pos3 | pre3 | pos4 | pre4 | pos5 | pre5 | pos6 | pre6 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ENSG00000000003.16 TSPAN6 chrX 100634689 100641991 | 10 | 69 | 119 | 3 | 44 | 48 | 91 | 126 | 61 | 83 | 73 | 25 |
| ENSG00000000005.6 TNMD chrX 100582936 100587019 | 12 | 37 | 67 | 30 | 17 | 56 | 69 | 84 | 62 | 65 | 43 | 24 |
| ENSG00000000005.6 TNMD chrX 100591624 100595624 | 39 | 32 | 40 | 54 | 30 | 46 | 36 | 34 | 34 | 42 | 27 | 18 |
| ENSG00000000419.14 DPM1 chr20 50940869 50951472 | 131 | 160 | 80 | 70 | 50 | 160 | 138 | 222 | 118 | 85 | 100 | 76 |
| ENSG00000000419.14 DPM1 chr20 50955153 50961140 | 40 | 40 | 64 | 9 | 78 | 45 | 170 | 173 | 103 | 170 | 107 | 89 |
| ENSG00000000457.14 SCYL3 chr1 169891896 169896267 | 56 | 61 | 102 | 92 | 44 | 69 | 101 | 77 | 57 | 144 | 73 | 91 |
colSums(k36a)## pos1 pre1 pos2 pre2 pos3 pre3 pos4 pre4
## 7230447 7422440 6514355 6192606 6754850 7790480 8050081 10962406
## pos5 pre5 pos6 pre6
## 7199883 8439267 6882750 8062324rpm9 <- apply(k9a, 2 , function(x) { x / sum(x) } ) * 1000000
rpm9 <- rpm9[rowMeans(rpm9) > 1,]
rpm36 <- apply(k36a, 2 , function(x) { x / sum(x) } ) * 1000000
rpm36 <- rpm36[rowMeans(rpm36) > 1,]MDS Plot
mds9 <- cmdscale(dist(t(k9a)))
cols=rep(c("pink","lightblue"),6)
XMIN=min(mds9[,1])*1.1
XMAX=max(mds9[,1])*1.1
plot(mds9, xlab="Coordinate 1", ylab="Coordinate 2", cex=2,col=cols,pch=19,main="H3K9ac",xlim=c(XMIN,XMAX),bty="n")
text(mds9, labels=colnames(k9a) ,cex=1.1)
pdf("k9a_mds.pdf",width=5,height=5)
plot(mds9, xlab="Coordinate 1", ylab="Coordinate 2", cex=2,col=cols,pch=19,main="H3K9ac",xlim=c(XMIN,XMAX),bty="n")
text(mds9, labels=colnames(k9a) ,cex=1.1)
dev.off()## png
## 2mds36 <- cmdscale(dist(t(k36a)))
cols=rep(c("pink","lightblue"),6)
XMIN=min(mds36[,1])*1.1
XMAX=max(mds36[,1])*1.1
plot(mds36, xlab="Coordinate 1", ylab="Coordinate 2", cex=2,col=cols,pch=19,main="H3K36ac",xlim=c(XMIN,XMAX),bty="n")
text(mds36, labels=colnames(k9a) ,cex=1.1)
pdf("k36a_mds.pdf",width=5,height=5)
plot(mds36, xlab="Coordinate 1", ylab="Coordinate 2", cex=2,col=cols,pch=19,main="H3K36ac",xlim=c(XMIN,XMAX),bty="n")
text(mds36, labels=colnames(k9a) ,cex=1.1)
dev.off()## png
## 2# include input sample
allac <- cbind(k9a,k36a)
mycolnames <- colnames(allac)
colnames(allac) <- paste(c(rep("k9",12),rep("k36",12)),mycolnames)
allac <- cbind(allac,input)
mdsa <- cmdscale(dist(t(allac)))
cols=c(rep(c("pink","lightblue"),12),rep("gray",4))
XMIN=min(mdsa[,1])*1.1
XMAX=max(mdsa[,1])*1.1
plot(mdsa, xlab="Coordinate 1", ylab="Coordinate 2", cex=2,col=cols,pch=19,main="ChIP",xlim=c(XMIN,XMAX),bty="n")
text(mdsa, labels=colnames(allac) ,cex=1.1)
# Remove input
allac <- cbind(k9a,k36a)
mycolnames <- colnames(allac)
colnames(allac) <- paste(c(rep("k9",12),rep("k36",12)),mycolnames)
mdsa <- cmdscale(dist(t(allac)))
cols=rep(c("pink","lightblue"),12)
XMIN=min(mdsa[,1])*1.1
XMAX=max(mdsa[,1])*1.1
plot(mdsa, xlab="Coordinate 1", ylab="Coordinate 2", cex=2,col=cols,pch=19,main="ChIP",xlim=c(XMIN,XMAX),bty="n")
text(mdsa, labels=colnames(allac) ,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(k9a))
ss$pos <- factor(grepl("pos",ss[,1]))
ss$participant <- c("1","1","2","2","3","3","4","4","5","5","6","6")
dim(k9a)## [1] 95874 12k9a <- k9a[rowMeans(k9a)>=10,]
dim(k9a)## [1] 92879 12dds <- DESeqDataSetFromMatrix(countData = k9a , colData = ss, design = ~ participant + pos )## 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))
dk9 <- as.data.frame(zz[order(zz$pvalue),])
write.table(dk9,file="k9a_deseq.tsv",quote=F,sep="\t")
head(subset(dk9,log2FoldChange > 0),20)[,c(1,2,6)] |>
kbl(caption = "Top increased genes in paired analysis for H3K9ac") |>
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | padj | |
|---|---|---|---|
| ENSG00000036672.16 USP2 chr11 119374356 119383711 | 524.6631 | 1.487334 | 0 |
| ENSG00000109819.9 PPARGC1A chr4 23887945 23892077 | 381.7201 | 1.567396 | 0 |
| ENSG00000173575.24 CHD2 chr15 92898189 92906908 | 484.2652 | 1.437051 | 0 |
| ENSG00000273338.1 ENSG00000273338 chr1 78002554 78006554 | 251.1536 | 1.849550 | 0 |
| ENSG00000101608.13 MYL12A chr18 3250218 3257524 | 268.4371 | 1.487999 | 0 |
| ENSG00000162616.9 DNAJB4 chr1 78002935 78007111 | 280.7745 | 1.755524 | 0 |
| ENSG00000100628.12 ASB2 chr14 93972083 93978739 | 421.3027 | 1.643101 | 0 |
| ENSG00000112183.15 RBM24 chr6 17287855 17291969 | 226.5850 | 1.790212 | 0 |
| ENSG00000150991.16 UBC chr12 124909716 124919368 | 485.9735 | 1.608265 | 0 |
| ENSG00000255958.1 GABARAPL1-AS1 chr12 10212761 10216761 | 294.3751 | 1.404936 | 0 |
| ENSG00000168386.19 FILIP1L chr3 100112264 100116513 | 200.6929 | 1.549050 | 0 |
| ENSG00000122643.23 NT5C3A chr7 33038900 33042910 | 219.3906 | 1.734015 | 0 |
| ENSG00000185437.16 SH3BGR chr21 39449970 39454123 | 218.1856 | 1.542191 | 0 |
| ENSG00000284324.1 MIR3175 chr15 92902399 92906399 | 257.2706 | 1.543437 | 0 |
| ENSG00000251562.11 MALAT1 chr11 65495640 65506702 | 954.1195 | 1.506696 | 0 |
| ENSG00000170807.12 LMOD2 chr7 123653866 123657964 | 249.8334 | 1.472954 | 0 |
| ENSG00000119508.18 NR4A3 chr9 99819855 99823885 | 378.4525 | 1.544204 | 0 |
| ENSG00000113448.21 PDE4D chr5 59584375 59588499 | 309.6319 | 1.472067 | 0 |
| ENSG00000265345.1 MIR5188 chr12 124913547 124917547 | 325.7779 | 1.653016 | 0 |
| ENSG00000168575.11 SLC20A2 chr8 42499231 42504702 | 246.5869 | 1.358816 | 0 |
head(subset(dk9,log2FoldChange < 0),20)[,c(1,2,6)] |>
kbl(caption = "Top increased genes in paired analysis for H3K9ac") |>
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | padj | |
|---|---|---|---|
| ENSG00000224078.15 SNHG14 chr15 25204080 25210108 | 99.84421 | -1.3083026 | 0.00e+00 |
| ENSG00000199833.1 SNORD115-21 chr15 25206083 25210083 | 69.81318 | -1.3420231 | 2.00e-07 |
| ENSG00000258991.1 DUX4L19 chrY 11330329 11334329 | 251.75574 | -1.0878236 | 7.00e-07 |
| ENSG00000075043.21 KCNQ2 chr20 63406012 63415439 | 170.49865 | -1.0313197 | 4.20e-06 |
| ENSG00000200801.1 SNORD115-28 chr15 25220348 25224348 | 68.88609 | -1.1700731 | 4.90e-06 |
| ENSG00000259029.1 DUX4L18 chrY 11319557 11323557 | 238.54566 | -0.9068931 | 6.60e-06 |
| ENSG00000196557.14 CACNA1H chr16 1208608 1214925 | 100.74329 | -1.1008554 | 8.80e-06 |
| ENSG00000143494.16 VASH2 chr1 212959227 212963227 | 80.14866 | -1.1582411 | 9.30e-06 |
| ENSG00000199704.1 SNORD115-29 chr15 25221246 25225246 | 71.97597 | -1.1749336 | 1.19e-05 |
| ENSG00000101292.8 PROKR2 chr20 5314414 5318954 | 100.20647 | -1.0905800 | 1.21e-05 |
| ENSG00000202261.1 SNORD115-44 chr15 25248859 25252859 | 57.54762 | -1.3339151 | 1.35e-05 |
| ENSG00000201496.1 RN7SKP275 chr4 5456881 5460881 | 75.21991 | -1.1833411 | 1.74e-05 |
| ENSG00000169325.10 GUSBP12 chr7 57201597 57205597 | 63.37317 | -1.1216681 | 1.90e-05 |
| ENSG00000203805.11 PLPP4 chr10 120501579 120505579 | 85.86486 | -1.0818387 | 2.05e-05 |
| ENSG00000259282.5 SPATA8-AS1 chr15 96779654 96785312 | 122.58135 | -0.9462947 | 2.09e-05 |
| ENSG00000155066.16 PROM2 chr2 95272449 95281866 | 165.45028 | -0.8805711 | 2.11e-05 |
| ENSG00000201969.1 SNORD115-20 chr15 25204262 25208262 | 67.76409 | -1.2025575 | 2.13e-05 |
| ENSG00000260134.1 C1QL1P1 chr16 35572834 35576834 | 65.10995 | -1.2059305 | 2.25e-05 |
| ENSG00000265378.1 ENSG00000265378 chr18 78503165 78507165 | 60.87092 | -1.1806800 | 2.25e-05 |
| ENSG00000258567.1 DUX4L16 chrY 11304918 11308918 | 187.92800 | -0.9672420 | 3.14e-05 |
nrow(subset(dk9,padj<0.05 & log2FoldChange > 0))## [1] 12279nrow(subset(dk9,padj<0.05 & log2FoldChange < 0))## [1] 4395dk9[grep("SLC2A4",rownames(dk9)),c(1,2,5,6)] |>
kbl(caption = "GLUT4 TSS H3K9ac)") |>
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | pvalue | padj | |
|---|---|---|---|---|
| ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 | 267.4548 | 0.8139358 | 0.0001994 | 0.0030238 |
| ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 | 233.0312 | 0.0703409 | 0.7450239 | 0.8602644 |
dk9[grep("SLC2A",rownames(dk9)),c(1,2,5,6)] |>
kbl(caption = "GLUT transporters TSS H3K9ac)") |>
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | pvalue | padj | |
|---|---|---|---|---|
| ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 | 267.45485 | 0.8139358 | 0.0001994 | 0.0030238 |
| ENSG00000142583.18 SLC2A5 chr1 9067536 9075122 | 172.73600 | 0.6381898 | 0.0002636 | 0.0037477 |
| ENSG00000146411.6 SLC2A12 chr6 134050480 134054624 | 144.56711 | 0.7032469 | 0.0007043 | 0.0078798 |
| ENSG00000133460.20 SLC2A11 chr22 23854703 23860136 | 212.78075 | 0.5438167 | 0.0044675 | 0.0305393 |
| ENSG00000059804.17 SLC2A3 chr12 7934034 7938187 | 110.18393 | 0.5448942 | 0.0058367 | 0.0369365 |
| ENSG00000117394.24 SLC2A1 chr1 42927499 42933769 | 122.39498 | -0.4902829 | 0.0071624 | 0.0426871 |
| ENSG00000059804.17 SLC2A3 chr12 7917750 7925289 | 99.74252 | -0.4671657 | 0.0115570 | 0.0598594 |
| ENSG00000173262.12 SLC2A14 chr12 7830735 7834816 | 64.48003 | -0.5050898 | 0.0227083 | 0.0954830 |
| ENSG00000109667.12 SLC2A9 chr4 10027300 10031300 | 66.96057 | 0.5004813 | 0.0340207 | 0.1260044 |
| ENSG00000109667.12 SLC2A9 chr4 10017075 10023497 | 160.84233 | 0.4630155 | 0.0344426 | 0.1270098 |
| ENSG00000173262.12 SLC2A14 chr12 7861455 7865455 | 64.04563 | -0.4777673 | 0.0346338 | 0.1275215 |
| ENSG00000173262.12 SLC2A14 chr12 7870380 7875484 | 77.94010 | -0.4514423 | 0.0372788 | 0.1340464 |
| ENSG00000142583.18 SLC2A5 chr1 9037895 9041895 | 74.98537 | -0.4304027 | 0.0500987 | 0.1631869 |
| ENSG00000163581.14 SLC2A2 chr3 171024670 171028743 | 85.22052 | 0.3657697 | 0.0717270 | 0.2074463 |
| ENSG00000109667.12 SLC2A9 chr4 9885642 9889642 | 65.42534 | -0.4175300 | 0.0807162 | 0.2240336 |
| ENSG00000185031.6 SLC2A3P2 chr1 64984087 64988087 | 56.99060 | -0.3928979 | 0.1096012 | 0.2737422 |
| ENSG00000250413.1 SLC2A9-AS1 chr4 10004482 10008482 | 83.75998 | 0.3215309 | 0.1223829 | 0.2936262 |
| ENSG00000117394.24 SLC2A1 chr1 42950016 42954016 | 79.60871 | -0.3285558 | 0.1746223 | 0.3673722 |
| ENSG00000173262.12 SLC2A14 chr12 7889127 7893148 | 61.39839 | -0.2891091 | 0.2062892 | 0.4076756 |
| ENSG00000059804.17 SLC2A3 chr12 8017004 8021007 | 94.42340 | 0.2422359 | 0.2215203 | 0.4254615 |
| ENSG00000275768.1 SLC2AXP1 chr2 95195357 95199357 | 64.99273 | -0.2757756 | 0.2333236 | 0.4394464 |
| ENSG00000142583.18 SLC2A5 chr1 9086443 9090478 | 74.50608 | -0.2334344 | 0.2819794 | 0.4927464 |
| ENSG00000253861.1 SLC2A3P1 chr5 168551727 168555727 | 56.81002 | -0.2683819 | 0.3063217 | 0.5186367 |
| ENSG00000151229.13 SLC2A13 chr12 39869881 39873881 | 72.09885 | -0.2349659 | 0.3117303 | 0.5239755 |
| ENSG00000109667.12 SLC2A9 chr4 10052875 10056936 | 64.58814 | -0.1974429 | 0.3790624 | 0.5883114 |
| ENSG00000197496.6 SLC2A10 chr20 46707649 46711737 | 128.93541 | -0.1505041 | 0.4649025 | 0.6625902 |
| ENSG00000109667.12 SLC2A9 chr4 9824600 9828600 | 72.88661 | -0.1660901 | 0.5028677 | 0.6936753 |
| ENSG00000109667.12 SLC2A9 chr4 9906345 9910345 | 59.19697 | -0.1432424 | 0.5522555 | 0.7308169 |
| ENSG00000254954.1 SLC2A13P1 chr8 60139821 60143821 | 66.56583 | 0.1718459 | 0.5541247 | 0.7319301 |
| ENSG00000254088.1 SLC2A3P4 chr8 86501591 86505591 | 55.65896 | -0.1192612 | 0.6387493 | 0.7916942 |
| ENSG00000163581.14 SLC2A2 chr3 171005470 171009470 | 73.36293 | -0.1024520 | 0.6539062 | 0.8014747 |
| ENSG00000151229.13 SLC2A13 chr12 40103808 40108089 | 189.60020 | 0.0797274 | 0.6709739 | 0.8128539 |
| ENSG00000133460.20 SLC2A11 chr22 23873160 23884960 | 183.83832 | 0.0782214 | 0.6722285 | 0.8137517 |
| ENSG00000117394.24 SLC2A1 chr1 42941238 42945238 | 90.92739 | -0.0736564 | 0.7052536 | 0.8360689 |
| ENSG00000117394.24 SLC2A1 chr1 42956651 42960893 | 144.13653 | 0.0730659 | 0.7345002 | 0.8540266 |
| ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 | 233.03122 | 0.0703409 | 0.7450239 | 0.8602644 |
| ENSG00000227533.7 SLC2A1-DT chr1 42957049 42961100 | 142.57101 | 0.0716479 | 0.7469866 | 0.8612671 |
| ENSG00000109667.12 SLC2A9 chr4 10038210 10042388 | 76.99335 | 0.0642317 | 0.7685924 | 0.8745123 |
| ENSG00000197496.6 SLC2A10 chr20 46712467 46716467 | 74.35784 | -0.0593869 | 0.7756693 | 0.8788244 |
| ENSG00000151229.13 SLC2A13 chr12 39833760 39837760 | 52.07064 | -0.0600615 | 0.8185127 | 0.9039442 |
| ENSG00000160326.14 SLC2A6 chr9 133477059 133481127 | 106.87263 | -0.0355185 | 0.8641603 | 0.9300374 |
| ENSG00000136856.18 SLC2A8 chr9 127395138 127409404 | 413.73250 | 0.0220683 | 0.8966627 | 0.9474715 |
| ENSG00000163581.14 SLC2A2 chr3 171016623 171020623 | 70.30672 | 0.0193509 | 0.9368715 | 0.9688972 |
| ENSG00000197241.4 SLC2A7 chr1 9024345 9028423 | 62.10553 | 0.0108188 | 0.9616640 | 0.9814776 |
| ENSG00000059804.17 SLC2A3 chr12 7928035 7932055 | 54.87079 | 0.0104012 | 0.9642192 | 0.9827446 |
Smear plot
sig9 <- subset(dk9,padj<0.05)
NSIG=nrow(sig9)
NDOWN=nrow(subset(sig9,log2FoldChange<0))
NUP=nrow(subset(sig9,log2FoldChange>0))
NTOT=nrow(dk9)
HEADER=paste(NTOT,"TSSs detected;",NSIG,"@5%FDR;",NUP,"up",NDOWN,"down")
plot(log10(dk9$baseMean),dk9$log2FoldChange,cex=0.5,pch=19,col="darkgray",bty="n",
xlab="log10 basemean",ylab="log2 fold change",main="H3K9ac")
points(log10(sig9$baseMean),sig9$log2FoldChange,cex=0.5,pch=19,col="black")
mtext(HEADER)
pdf("k9a_smear.pdf",width=5,height=5)
plot(log10(dk9$baseMean),dk9$log2FoldChange,cex=0.5,pch=19,col="darkgray",btn="n",
xlab="log10 basemean",ylab="log2 fold change",main="H3K9ac")## Warning in plot.window(...): "btn" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "btn" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "btn" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "btn" is not a
## graphical parameter## Warning in box(...): "btn" is not a graphical parameter## Warning in title(...): "btn" is not a graphical parameterpoints(log10(sig9$baseMean),sig9$log2FoldChange,cex=0.5,pch=19,col="black")
mtext(HEADER)
dev.off()## png
## 2subset(dk9,log10(baseMean)>3)[,c(1,2,5,6)]|>
kbl(caption = "TSS with high basemeans") |>
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | pvalue | padj | |
|---|---|---|---|---|
| ENSG00000211459.2 MT-RNR1 chrM 1 2648 | 1531.180 | -0.6286146 | 0.1766221 | 0.3700204 |
| ENSG00000210049.1 MT-TF chrM 1 2577 | 1483.360 | -0.6233035 | 0.1826570 | 0.3777976 |
| ENSG00000228253.1 MT-ATP8 chrM 6366 10570 | 1782.739 | -0.6430670 | 0.1890921 | 0.3864007 |
| ENSG00000210156.1 MT-TK chrM 6295 10295 | 1656.142 | -0.6445407 | 0.1952980 | 0.3943193 |
| ENSG00000198899.2 MT-ATP6 chrM 6527 11205 | 2060.259 | -0.6400571 | 0.1958356 | 0.3950162 |
| ENSG00000210077.1 MT-TV chrM 1 3602 | 2166.439 | -0.6180571 | 0.1981148 | 0.3977412 |
| ENSG00000198938.2 MT-CO3 chrM 7207 11207 | 1874.265 | -0.6381460 | 0.2001226 | 0.4001461 |
| ENSG00000210082.2 MT-RNR2 chrM 1 3671 | 2213.631 | -0.6102104 | 0.2047901 | 0.4058425 |
| ENSG00000210196.2 MT-TP chrM 14023 16568 | 1993.267 | -0.6030131 | 0.2058250 | 0.4070699 |
| ENSG00000198712.1 MT-CO2 chrM 5586 10267 | 1806.188 | -0.6350773 | 0.2098888 | 0.4121003 |
| ENSG00000198695.2 MT-ND6 chrM 12673 16568 | 2924.634 | -0.6084969 | 0.2120392 | 0.4146772 |
| ENSG00000210194.1 MT-TE chrM 12742 16568 | 2885.734 | -0.6083971 | 0.2138021 | 0.4167588 |
| ENSG00000198727.2 MT-CYB chrM 12747 16568 | 2882.337 | -0.6081162 | 0.2138906 | 0.4168701 |
| ENSG00000210164.1 MT-TG chrM 7991 11991 | 2107.927 | -0.6270086 | 0.2157150 | 0.4189402 |
| ENSG00000210195.2 MT-TT chrM 13888 16568 | 2093.509 | -0.5854327 | 0.2171742 | 0.4207623 |
| ENSG00000198840.2 MT-ND3 chrM 8059 12059 | 2128.446 | -0.6197204 | 0.2213776 | 0.4253218 |
| ENSG00000210154.1 MT-TD chrM 5518 9518 | 1357.741 | -0.6282122 | 0.2223644 | 0.4264527 |
| ENSG00000210151.2 MT-TS1 chrM 5514 9514 | 1357.739 | -0.6273058 | 0.2233567 | 0.4277526 |
| ENSG00000210174.1 MT-TR chrM 8405 12405 | 2234.079 | -0.6066202 | 0.2321942 | 0.4381010 |
| ENSG00000209082.1 MT-TL1 chrM 1230 5230 | 1815.209 | -0.5913926 | 0.2365977 | 0.4432577 |
| ENSG00000198888.2 MT-ND1 chrM 1307 5307 | 1818.061 | -0.5875696 | 0.2380542 | 0.4450083 |
| ENSG00000198886.2 MT-ND4 chrM 8760 12760 | 2334.944 | -0.5902927 | 0.2465821 | 0.4545098 |
| ENSG00000210176.1 MT-TH chrM 10138 14138 | 2700.257 | -0.5865324 | 0.2490944 | 0.4570092 |
| ENSG00000210184.1 MT-TS2 chrM 10207 14207 | 2681.210 | -0.5846668 | 0.2499794 | 0.4580509 |
| ENSG00000212907.2 MT-ND4L chrM 8470 12764 | 2456.237 | -0.5761765 | 0.2535204 | 0.4622171 |
| ENSG00000210191.1 MT-TL2 chrM 10266 14266 | 2665.280 | -0.5803521 | 0.2542443 | 0.4629461 |
| ENSG00000210140.1 MT-TC chrM 3826 7826 | 1272.745 | -0.5695709 | 0.2677136 | 0.4779111 |
| ENSG00000198804.2 MT-CO1 chrM 3904 7904 | 1307.687 | -0.5677777 | 0.2686811 | 0.4789795 |
| ENSG00000210144.1 MT-TY chrM 3891 7891 | 1309.665 | -0.5665567 | 0.2691696 | 0.4795282 |
| ENSG00000210100.1 MT-TI chrM 2263 6263 | 1753.693 | -0.5664499 | 0.2795816 | 0.4902671 |
| ENSG00000210117.1 MT-TW chrM 3512 7512 | 1374.602 | -0.5587773 | 0.2806861 | 0.4914401 |
| ENSG00000210112.1 MT-TM chrM 2402 6402 | 1752.737 | -0.5674392 | 0.2822875 | 0.4930331 |
| ENSG00000210107.1 MT-TQ chrM 2400 6400 | 1752.365 | -0.5665396 | 0.2827170 | 0.4934226 |
| ENSG00000210135.1 MT-TN chrM 3729 7729 | 1300.006 | -0.5467415 | 0.2866848 | 0.4977334 |
| ENSG00000198763.3 MT-ND2 chrM 2470 6470 | 1746.583 | -0.5585899 | 0.2897799 | 0.5009393 |
| ENSG00000198786.2 MT-ND5 chrM 10337 16146 | 3964.289 | -0.5353781 | 0.2918433 | 0.5030967 |
| ENSG00000210127.1 MT-TA chrM 3655 7655 | 1315.246 | -0.5329101 | 0.3002276 | 0.5119965 |
Now for H3K36a
dim(k36a)## [1] 95874 12k36a <- k36a[rowMeans(k36a)>=10,]
dim(k36a)## [1] 92792 12dds <- DESeqDataSetFromMatrix(countData = k36a , colData = ss, design = ~ participant + pos )## 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## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.## final dispersion estimates## fitting model and testingz <- results(res)
vsd <- vst(dds, blind=FALSE)## -- note: fitType='parametric', but the dispersion trend was not well captured by the
## function: y = a/x + b, and a local regression fit was automatically substituted.
## specify fitType='local' or 'mean' to avoid this message next time.zz <- cbind(as.data.frame(z),assay(vsd))
dk36 <- as.data.frame(zz[order(zz$pvalue),])
write.table(dk36,file="k36a_deseq.tsv",quote=F,sep="\t")
head(subset(dk36,log2FoldChange > 0),20)[,c(1,2,6)] |>
kbl(caption = "Top increased genes in paired analysis for H3K36ac") |>
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | padj | |
|---|---|---|---|
| ENSG00000096060.15 FKBP5 chr6 35726583 35730583 | 226.4563 | 1.734588 | 0 |
| ENSG00000072422.17 RHOBTB1 chr10 60942185 60946185 | 159.8681 | 1.924122 | 0 |
| ENSG00000235831.8 BHLHE40-AS1 chr3 4977961 4982414 | 249.6403 | 1.801638 | 0 |
| ENSG00000134107.5 BHLHE40 chr3 4977437 4982232 | 273.0280 | 1.775006 | 0 |
| ENSG00000177354.12 C10orf71 chr10 49297170 49301274 | 182.2851 | 1.877299 | 0 |
| ENSG00000289989.1 ENSG00000289989 chr10 60942378 60946378 | 162.3501 | 1.807843 | 0 |
| ENSG00000150991.16 UBC chr12 124909716 124919368 | 277.3492 | 1.405234 | 0 |
| ENSG00000290057.1 ENSG00000290057 chr11 65419448 65423448 | 199.8796 | 1.938683 | 0 |
| ENSG00000154803.13 FLCN chr17 17235139 17239188 | 122.8297 | 1.896609 | 0 |
| ENSG00000264187.1 ENSG00000264187 chr17 17235185 17239185 | 122.8297 | 1.896609 | 0 |
| ENSG00000166272.18 WBP1L chr10 102774249 102778391 | 156.6354 | 2.207488 | 0 |
| ENSG00000259164.1 ENSG00000259164 chr14 77096311 77100311 | 140.4952 | 2.185755 | 0 |
| ENSG00000266498.1 ENSG00000266498 chr17 17234118 17238118 | 122.5840 | 1.729711 | 0 |
| ENSG00000236208.3 C10orf71-AS1 chr10 49296775 49301018 | 182.3681 | 1.869396 | 0 |
| ENSG00000198894.8 CIPC chr14 77096126 77101176 | 166.0308 | 2.116770 | 0 |
| ENSG00000227630.4 LINC01132 chr1 234722042 234727398 | 152.0404 | 1.680520 | 0 |
| ENSG00000167716.19 WDR81 chr17 1714523 1718601 | 175.3965 | 1.735130 | 0 |
| ENSG00000248538.10 PPP1R3B-DT chr8 9149695 9157290 | 217.0624 | 1.594199 | 0 |
| ENSG00000054523.20 KIF1B chr1 10208570 10213616 | 148.4638 | 1.621419 | 0 |
| ENSG00000237461.3 ENSG00000237461 chr9 99817889 99821921 | 163.7983 | 1.882090 | 0 |
head(subset(dk36,log2FoldChange < 0),20)[,c(1,2,6)] |>
kbl(caption = "Top increased genes in paired analysis for H3K36ac") |>
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | padj | |
|---|---|---|---|
| ENSG00000267073.1 ENSG00000267073 chr19 1746056 1750160 | 86.43211 | -2.415936 | 0 |
| ENSG00000176884.17 GRIN1 chr9 137160706 137169923 | 215.89789 | -2.151658 | 0 |
| ENSG00000182871.16 COL18A1 chr21 45503358 45508019 | 107.17483 | -1.812416 | 0 |
| ENSG00000203900.2 ENSG00000203900 chr20 63357988 63361988 | 80.28738 | -2.261360 | 0 |
| ENSG00000174407.15 MIR1-1HG chr20 62568442 62572442 | 76.48334 | -2.479432 | 0 |
| ENSG00000261678.3 SCRT1 chr8 144334169 144338482 | 74.54850 | -3.002328 | 0 |
| ENSG00000141577.14 CEP131 chr17 81190376 81199019 | 216.21803 | -1.328939 | 0 |
| ENSG00000125503.13 PPP1R12C chr19 55090598 55099064 | 117.73322 | -2.029696 | 0 |
| ENSG00000261793.1 ENSG00000261793 chr9 137169984 137173984 | 71.86733 | -2.223339 | 0 |
| ENSG00000007516.15 BAIAP3 chr16 1339822 1349298 | 196.29268 | -1.614631 | 0 |
| ENSG00000188064.10 WNT7B chr22 45970237 45979162 | 187.13948 | -1.481744 | 0 |
| ENSG00000162004.19 CCDC78 chr16 722635 728954 | 141.26145 | -1.860519 | 0 |
| ENSG00000203485.14 INF2 chr14 104697368 104717725 | 428.49144 | -1.352265 | 0 |
| ENSG00000177106.17 EPS8L2 chr11 718606 728954 | 217.68061 | -1.699419 | 0 |
| ENSG00000099821.14 POLRMT chr19 615997 621226 | 63.51642 | -2.389155 | 0 |
| ENSG00000165716.11 DIPK1B chr9 136718885 136723984 | 121.69781 | -1.877862 | 0 |
| ENSG00000105298.14 CACTIN chr19 3609925 3623028 | 240.82412 | -1.339066 | 0 |
| ENSG00000127561.16 SYNGR3 chr16 1987660 1994463 | 134.80975 | -1.396684 | 0 |
| ENSG00000259803.8 SLC22A31 chr16 89196840 89203651 | 125.36478 | -1.892523 | 0 |
| ENSG00000174938.15 SEZ6L2 chr16 29895928 29901550 | 108.38265 | -1.926817 | 0 |
nrow(subset(dk36,padj<0.05 & log2FoldChange > 0))## [1] 2885nrow(subset(dk36,padj<0.05 & log2FoldChange < 0))## [1] 4230dk36[grep("SLC2A4",rownames(dk36)),c(1,2,5,6)] |>
kbl(caption = "GLUT4 TSS H3K36ac)") |>
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | pvalue | padj | |
|---|---|---|---|---|
| ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 | 186.6098 | 1.3001087 | 0.0000000 | 0.0000003 |
| ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 | 135.0512 | -0.4546132 | 0.0392601 | 0.2045954 |
dk36[grep("SLC2A",rownames(dk36)),c(1,2,5,6)] |>
kbl(caption = "GLUT transporters TSS H3K36ac)") |>
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | pvalue | padj | |
|---|---|---|---|---|
| ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 | 186.60980 | 1.3001087 | 0.0000000 | 0.0000003 |
| ENSG00000146411.6 SLC2A12 chr6 134050480 134054624 | 75.05498 | 1.1507533 | 0.0001220 | 0.0047234 |
| ENSG00000142583.18 SLC2A5 chr1 9037895 9041895 | 85.52562 | -1.1976609 | 0.0003074 | 0.0092344 |
| ENSG00000254954.1 SLC2A13P1 chr8 60139821 60143821 | 62.95483 | 1.2452427 | 0.0003587 | 0.0102570 |
| ENSG00000227533.7 SLC2A1-DT chr1 42957049 42961100 | 85.88964 | -1.0391178 | 0.0007386 | 0.0168361 |
| ENSG00000117394.24 SLC2A1 chr1 42956651 42960893 | 88.03535 | -1.0263314 | 0.0008009 | 0.0177798 |
| ENSG00000173262.12 SLC2A14 chr12 7870380 7875484 | 94.36016 | -0.8382695 | 0.0043144 | 0.0539400 |
| ENSG00000275768.1 SLC2AXP1 chr2 95195357 95199357 | 61.71953 | 1.1725857 | 0.0067826 | 0.0720183 |
| ENSG00000136856.18 SLC2A8 chr9 127395138 127409404 | 263.42910 | -0.6147938 | 0.0070274 | 0.0737636 |
| ENSG00000142583.18 SLC2A5 chr1 9067536 9075122 | 186.21955 | 0.5119162 | 0.0118638 | 0.1017092 |
| ENSG00000163581.14 SLC2A2 chr3 171005470 171009470 | 83.61251 | 0.5699882 | 0.0250311 | 0.1585450 |
| ENSG00000160326.14 SLC2A6 chr9 133477059 133481127 | 63.26441 | -0.8467236 | 0.0338489 | 0.1878643 |
| ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 | 135.05116 | -0.4546132 | 0.0392601 | 0.2045954 |
| ENSG00000151229.13 SLC2A13 chr12 40103808 40108089 | 72.00256 | -0.6898400 | 0.0459551 | 0.2239637 |
| ENSG00000185031.6 SLC2A3P2 chr1 64984087 64988087 | 70.08017 | 0.7549150 | 0.0504972 | 0.2362238 |
| ENSG00000173262.12 SLC2A14 chr12 7861455 7865455 | 79.27207 | -0.4934816 | 0.0603866 | 0.2609641 |
| ENSG00000253861.1 SLC2A3P1 chr5 168551727 168555727 | 58.64129 | 0.5568002 | 0.0827752 | 0.3112368 |
| ENSG00000250413.1 SLC2A9-AS1 chr4 10004482 10008482 | 74.65561 | 0.4746595 | 0.0988203 | 0.3430500 |
| ENSG00000163581.14 SLC2A2 chr3 171016623 171020623 | 76.09859 | 0.4292654 | 0.1462781 | 0.4206863 |
| ENSG00000133460.20 SLC2A11 chr22 23873160 23884960 | 188.19388 | -0.3491312 | 0.1677292 | 0.4511284 |
| ENSG00000117394.24 SLC2A1 chr1 42941238 42945238 | 85.58128 | -0.3984973 | 0.2123230 | 0.5087292 |
| ENSG00000151229.13 SLC2A13 chr12 39869881 39873881 | 56.53905 | 0.5333175 | 0.2136755 | 0.5103357 |
| ENSG00000197241.4 SLC2A7 chr1 9024345 9028423 | 70.94595 | -0.4711149 | 0.2203441 | 0.5185961 |
| ENSG00000133460.20 SLC2A11 chr22 23854703 23860136 | 57.86951 | -0.3740701 | 0.2796223 | 0.5794585 |
| ENSG00000109667.12 SLC2A9 chr4 9824600 9828600 | 59.67080 | 0.3267420 | 0.3190199 | 0.6166880 |
| ENSG00000197496.6 SLC2A10 chr20 46707649 46711737 | 84.94141 | -0.2230179 | 0.3984021 | 0.6838454 |
| ENSG00000109667.12 SLC2A9 chr4 10017075 10023497 | 135.99731 | 0.2143309 | 0.4142911 | 0.6965651 |
| ENSG00000109667.12 SLC2A9 chr4 10027300 10031300 | 91.57337 | -0.3141933 | 0.4243097 | 0.7043559 |
| ENSG00000059804.17 SLC2A3 chr12 7917750 7925289 | 91.26654 | -0.2104775 | 0.4397228 | 0.7156218 |
| ENSG00000109667.12 SLC2A9 chr4 9906345 9910345 | 65.75939 | -0.3412816 | 0.4455980 | 0.7198539 |
| ENSG00000109667.12 SLC2A9 chr4 10038210 10042388 | 60.53255 | -0.3370922 | 0.4691768 | 0.7362232 |
| ENSG00000163581.14 SLC2A2 chr3 171024670 171028743 | 63.91328 | -0.2226109 | 0.5262352 | 0.7746064 |
| ENSG00000117394.24 SLC2A1 chr1 42950016 42954016 | 67.78232 | -0.1772228 | 0.5342280 | 0.7798849 |
| ENSG00000059804.17 SLC2A3 chr12 8017004 8021007 | 53.73041 | 0.1806048 | 0.6161042 | 0.8292756 |
| ENSG00000059804.17 SLC2A3 chr12 7934034 7938187 | 66.10362 | -0.1025274 | 0.7183913 | 0.8831709 |
| ENSG00000109667.12 SLC2A9 chr4 10052875 10056936 | 63.69302 | 0.2513374 | 0.7589396 | 0.9034469 |
| ENSG00000059804.17 SLC2A3 chr12 7928035 7932055 | 58.09286 | -0.1137070 | 0.7710656 | 0.9088642 |
| ENSG00000254088.1 SLC2A3P4 chr8 86501591 86505591 | 53.71431 | 0.0784831 | 0.8475069 | 0.9421950 |
| ENSG00000109667.12 SLC2A9 chr4 9885642 9889642 | 83.21588 | -0.0398711 | 0.8975956 | 0.9616858 |
| ENSG00000151229.13 SLC2A13 chr12 39833760 39837760 | 34.50660 | 0.0617107 | 0.9119888 | 0.9676323 |
| ENSG00000197496.6 SLC2A10 chr20 46712467 46716467 | 97.04291 | 0.0296196 | 0.9406155 | 0.9780218 |
| ENSG00000142583.18 SLC2A5 chr1 9086443 9090478 | 59.87187 | 0.0314576 | 0.9431996 | 0.9791046 |
| ENSG00000117394.24 SLC2A1 chr1 42927499 42933769 | 123.94251 | -0.0045763 | 0.9827172 | 0.9940423 |
| ENSG00000173262.12 SLC2A14 chr12 7889127 7893148 | 57.22430 | -0.0062789 | 0.9848379 | 0.9948700 |
| ENSG00000173262.12 SLC2A14 chr12 7830735 7834816 | 69.97344 | -0.0026080 | 0.9929197 | 0.9976072 |
Smear plot
sig36 <- subset(dk36,padj<0.05)
NSIG=nrow(sig36)
NDOWN=nrow(subset(sig36,log2FoldChange<0))
NUP=nrow(subset(sig36,log2FoldChange>0))
NTOT=nrow(dk36)
HEADER=paste(NTOT,"TSSs detected;",NSIG,"@5%FDR;",NUP,"up",NDOWN,"down")
plot(log10(dk36$baseMean),dk36$log2FoldChange,cex=0.5,pch=19,col="darkgray",
xlab="log10 basemean",ylab="log2 fold change",main="H3K36ac")
points(log10(sig36$baseMean),sig36$log2FoldChange,cex=0.5,pch=19,col="black")
mtext(HEADER)
pdf("k36a_smear.pdf",width=5,height=5)
plot(log10(dk36$baseMean),dk36$log2FoldChange,cex=0.5,pch=19,col="darkgray",
xlab="log10 basemean",ylab="log2 fold change",main="H3K36ac")
points(log10(sig36$baseMean),sig36$log2FoldChange,cex=0.5,pch=19,col="black")
mtext(HEADER)
dev.off()## png
## 2subset(dk36,log10(baseMean)>3)[,c(1,2,5,6)] |>
kbl(caption = "TSS with high basemeans") |>
kable_paper("hover", full_width = F)| baseMean | log2FoldChange | pvalue | padj | |
|---|---|---|---|---|
| ENSG00000210127.1 MT-TA chrM 3655 7655 | 2135.902 | -2.343822 | 0.0000000 | 0.0000024 |
| ENSG00000210135.1 MT-TN chrM 3729 7729 | 2076.698 | -2.316781 | 0.0000000 | 0.0000029 |
| ENSG00000210140.1 MT-TC chrM 3826 7826 | 2006.645 | -2.339622 | 0.0000000 | 0.0000037 |
| ENSG00000210144.1 MT-TY chrM 3891 7891 | 2047.510 | -2.330088 | 0.0000000 | 0.0000040 |
| ENSG00000198804.2 MT-CO1 chrM 3904 7904 | 2044.241 | -2.307781 | 0.0000000 | 0.0000049 |
| ENSG00000210151.2 MT-TS1 chrM 5514 9514 | 2129.082 | -2.251115 | 0.0000000 | 0.0000061 |
| ENSG00000210117.1 MT-TW chrM 3512 7512 | 2292.611 | -2.308491 | 0.0000000 | 0.0000064 |
| ENSG00000210154.1 MT-TD chrM 5518 9518 | 2122.247 | -2.244491 | 0.0000000 | 0.0000066 |
| ENSG00000210049.1 MT-TF chrM 1 2577 | 6434.538 | -2.926481 | 0.0000007 | 0.0001111 |
| ENSG00000198712.1 MT-CO2 chrM 5586 10267 | 2781.168 | -2.171718 | 0.0000009 | 0.0001331 |
| ENSG00000211459.2 MT-RNR1 chrM 1 2648 | 6516.761 | -2.913611 | 0.0000009 | 0.0001394 |
| ENSG00000210156.1 MT-TK chrM 6295 10295 | 2553.785 | -2.150389 | 0.0000010 | 0.0001497 |
| ENSG00000210196.2 MT-TP chrM 14023 16568 | 8364.401 | -3.061226 | 0.0000013 | 0.0001838 |
| ENSG00000228253.1 MT-ATP8 chrM 6366 10570 | 2750.381 | -2.129767 | 0.0000016 | 0.0002134 |
| ENSG00000210195.2 MT-TT chrM 13888 16568 | 8548.123 | -3.018129 | 0.0000021 | 0.0002550 |
| ENSG00000210100.1 MT-TI chrM 2263 6263 | 3030.833 | -2.192392 | 0.0000025 | 0.0002916 |
| ENSG00000210107.1 MT-TQ chrM 2400 6400 | 3065.603 | -2.224794 | 0.0000025 | 0.0002956 |
| ENSG00000210112.1 MT-TM chrM 2402 6402 | 3065.732 | -2.224896 | 0.0000026 | 0.0002974 |
| ENSG00000198763.3 MT-ND2 chrM 2470 6470 | 3072.559 | -2.221097 | 0.0000032 | 0.0003525 |
| ENSG00000209082.1 MT-TL1 chrM 1230 5230 | 3062.922 | -2.189319 | 0.0000045 | 0.0004480 |
| ENSG00000198888.2 MT-ND1 chrM 1307 5307 | 3082.437 | -2.190628 | 0.0000047 | 0.0004673 |
| ENSG00000198938.2 MT-CO3 chrM 7207 11207 | 2950.215 | -2.035901 | 0.0000068 | 0.0006144 |
| ENSG00000198899.2 MT-ATP6 chrM 6527 11205 | 3197.516 | -2.065793 | 0.0000069 | 0.0006159 |
| ENSG00000210077.1 MT-TV chrM 1 3602 | 7720.904 | -2.764251 | 0.0000102 | 0.0008106 |
| ENSG00000210082.2 MT-RNR2 chrM 1 3671 | 7824.990 | -2.764046 | 0.0000110 | 0.0008579 |
| ENSG00000210194.1 MT-TE chrM 12742 16568 | 9858.920 | -2.914198 | 0.0000139 | 0.0010039 |
| ENSG00000198727.2 MT-CYB chrM 12747 16568 | 9852.251 | -2.912915 | 0.0000142 | 0.0010159 |
| ENSG00000198695.2 MT-ND6 chrM 12673 16568 | 9927.329 | -2.896127 | 0.0000163 | 0.0011187 |
| ENSG00000210174.1 MT-TR chrM 8405 12405 | 3621.956 | -2.033733 | 0.0000260 | 0.0015844 |
| ENSG00000198840.2 MT-ND3 chrM 8059 12059 | 3430.192 | -2.007427 | 0.0000352 | 0.0019477 |
| ENSG00000210164.1 MT-TG chrM 7991 11991 | 3393.323 | -1.996778 | 0.0000368 | 0.0020155 |
| ENSG00000198886.2 MT-ND4 chrM 8760 12760 | 3744.990 | -2.015837 | 0.0000462 | 0.0023618 |
| ENSG00000212907.2 MT-ND4L chrM 8470 12764 | 3966.438 | -2.008303 | 0.0000598 | 0.0028498 |
| ENSG00000210191.1 MT-TL2 chrM 10266 14266 | 4327.533 | -2.088604 | 0.0000914 | 0.0038490 |
| ENSG00000210184.1 MT-TS2 chrM 10207 14207 | 4351.605 | -2.091032 | 0.0000931 | 0.0038962 |
| ENSG00000210176.1 MT-TH chrM 10138 14138 | 4390.066 | -2.080087 | 0.0001088 | 0.0043462 |
| ENSG00000198786.2 MT-ND5 chrM 10337 16146 | 7388.176 | -2.265596 | 0.0001958 | 0.0066978 |
Mitch
H3K9ac
pw <- gmt_import("ref/c5.go.v2024.1.Hs.symbols.gmt")
names(pw) <- gsub("_", " ", names(pw))
gt2 <- as.data.frame(rownames(k9a))
gt2$gene <- sapply(strsplit(gt2[,1]," "),"[[",2)
m9 <- mitch_import(x=dk9,DEtype="deseq2",geneTable=gt2)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 92879## Note: no. genes in output = 59180## Note: estimated proportion of input genes in output = 0.637head(m9)## x
## 5_8S_rRNA -3.250775
## 5S_rRNA -1.211059
## 7SK 1.031242
## A1BG -2.304353
## A1BG-AS1 -1.932876
## A1CF -1.500519mres9 <- mitch_calc(x=m9,genesets=pw,minsetsize=5,cores=8,priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(subset(mres9$enrichment_result,p.adjustANOVA<0.05 & s.dist>0),50) |>
kbl(caption = "Top upregulated GO terms based on effect size (FDR<0.05) for H3K9ac") |>
kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 8557 | GOCC TRANSCRIPTION FACTOR AP 1 COMPLEX | 5 | 0.0004106 | 0.9122738 | 0.0014200 |
| 9159 | GOMF FATZ BINDING | 5 | 0.0004607 | 0.9043853 | 0.0015704 |
| 2846 | GOBP NEGATIVE REGULATION OF CALCIUM ION EXPORT ACROSS PLASMA MEMBRANE | 5 | 0.0004748 | 0.9023033 | 0.0016149 |
| 8449 | GOCC SARCOPLASMIC RETICULUM LUMEN | 9 | 0.0000040 | 0.8871707 | 0.0000206 |
| 8108 | GOCC MCRD MEDIATED MRNA STABILITY COMPLEX | 5 | 0.0006419 | 0.8813080 | 0.0021177 |
| 7842 | GOCC CYTOPLASMIC SIDE OF ROUGH ENDOPLASMIC RETICULUM MEMBRANE | 5 | 0.0006770 | 0.8775564 | 0.0022147 |
| 4996 | GOBP PROTEIN INSERTION INTO MEMBRANE FROM INNER SIDE | 5 | 0.0007410 | 0.8711618 | 0.0023974 |
| 8095 | GOCC LSM1 7 PAT1 COMPLEX | 5 | 0.0007665 | 0.8687554 | 0.0024684 |
| 4998 | GOBP PROTEIN INSERTION INTO MITOCHONDRIAL MEMBRANE | 7 | 0.0000717 | 0.8664883 | 0.0002923 |
| 8596 | GOCC U6 SNRNP | 11 | 0.0000009 | 0.8540851 | 0.0000052 |
| 8096 | GOCC LSM2 8 COMPLEX | 7 | 0.0000928 | 0.8529783 | 0.0003716 |
| 10170 | GOMF SMALL RIBOSOMAL SUBUNIT RRNA BINDING | 8 | 0.0000295 | 0.8527217 | 0.0001287 |
| 8147 | GOCC MITOCHONDRIAL INTERMEMBRANE SPACE PROTEIN TRANSPORTER COMPLEX | 6 | 0.0003125 | 0.8496581 | 0.0011136 |
| 8813 | GOMF ATP DEPENDENT PROTEIN BINDING | 6 | 0.0003128 | 0.8496017 | 0.0011142 |
| 8351 | GOCC PREFOLDIN COMPLEX | 7 | 0.0000992 | 0.8494391 | 0.0003945 |
| 3648 | GOBP NUCLEAR TRANSCRIBED MRNA CATABOLIC PROCESS DEADENYLATION INDEPENDENT DECAY | 7 | 0.0001124 | 0.8427806 | 0.0004420 |
| 5385 | GOBP REGULATION OF CALCIUM ION EXPORT ACROSS PLASMA MEMBRANE | 6 | 0.0003636 | 0.8403353 | 0.0012753 |
| 5597 | GOBP REGULATION OF ELECTRON TRANSFER ACTIVITY | 5 | 0.0011375 | 0.8402231 | 0.0035052 |
| 10199 | GOMF SPLICING FACTOR BINDING | 6 | 0.0003690 | 0.8394284 | 0.0012911 |
| 7887 | GOCC EKC KEOPS COMPLEX | 5 | 0.0011559 | 0.8390469 | 0.0035566 |
| 2997 | GOBP NEGATIVE REGULATION OF EPITHELIAL CELL PROLIFERATION INVOLVED IN PROSTATE GLAND DEVELOPMENT | 6 | 0.0003951 | 0.8351866 | 0.0013752 |
| 10363 | GOMF U6 SNRNA BINDING | 16 | 0.0000000 | 0.8317537 | 0.0000001 |
| 7794 | GOCC COMMITMENT COMPLEX | 5 | 0.0013364 | 0.8283329 | 0.0040498 |
| 8404 | GOCC RAGULATOR COMPLEX | 6 | 0.0004567 | 0.8261342 | 0.0015590 |
| 7713 | GOCC CCAAT BINDING FACTOR COMPLEX | 5 | 0.0014027 | 0.8247300 | 0.0042311 |
| 7359 | GOBP TRNA THREONYLCARBAMOYLADENOSINE MODIFICATION | 5 | 0.0014594 | 0.8217693 | 0.0043820 |
| 1588 | GOBP FORMATION OF QUADRUPLE SL U4 U5 U6 SNRNP | 10 | 0.0000069 | 0.8212439 | 0.0000334 |
| 7580 | GOBP WYBUTOSINE METABOLIC PROCESS | 6 | 0.0005601 | 0.8132400 | 0.0018709 |
| 8582 | GOCC TRNA METHYLTRANSFERASE COMPLEX | 7 | 0.0001998 | 0.8116419 | 0.0007441 |
| 10370 | GOMF UBIQUITIN LIGASE INHIBITOR ACTIVITY | 9 | 0.0000326 | 0.7995227 | 0.0001411 |
| 9765 | GOMF OXIDATIVE PHOSPHORYLATION UNCOUPLER ACTIVITY | 5 | 0.0020127 | 0.7974280 | 0.0058076 |
| 9975 | GOMF PROTEASOME ACTIVATING ACTIVITY | 6 | 0.0007258 | 0.7966055 | 0.0023542 |
| 7926 | GOCC EUKARYOTIC TRANSLATION INITIATION FACTOR 2B COMPLEX | 5 | 0.0020439 | 0.7962450 | 0.0058863 |
| 9943 | GOMF PORIN ACTIVITY | 6 | 0.0007722 | 0.7925835 | 0.0024838 |
| 2569 | GOBP MITOCHONDRIAL PROTON TRANSPORTING ATP SYNTHASE COMPLEX ASSEMBLY | 8 | 0.0001038 | 0.7923258 | 0.0004108 |
| 7911 | GOCC ENDOPLASMIC RETICULUM TUBULAR NETWORK MEMBRANE | 5 | 0.0022035 | 0.7904520 | 0.0062732 |
| 2580 | GOBP MITOCHONDRIAL TRANSLATIONAL ELONGATION | 5 | 0.0023117 | 0.7867410 | 0.0065423 |
| 6553 | GOBP RESPIRATORY CHAIN COMPLEX III ASSEMBLY | 12 | 0.0000024 | 0.7861288 | 0.0000128 |
| 1214 | GOBP DORSAL AORTA MORPHOGENESIS | 8 | 0.0001266 | 0.7824224 | 0.0004931 |
| 8802 | GOMF ATPASE INHIBITOR ACTIVITY | 5 | 0.0024694 | 0.7816037 | 0.0069158 |
| 890 | GOBP CILIARY NEUROTROPHIC FACTOR MEDIATED SIGNALING PATHWAY | 5 | 0.0024974 | 0.7807250 | 0.0069877 |
| 5052 | GOBP PROTEIN LOCALIZATION TO NUCLEAR BODY | 12 | 0.0000031 | 0.7768304 | 0.0000165 |
| 5114 | GOBP PROTEIN TO MEMBRANE DOCKING | 5 | 0.0026270 | 0.7767638 | 0.0073014 |
| 7846 | GOCC CYTOSOLIC LARGE RIBOSOMAL SUBUNIT | 57 | 0.0000000 | 0.7765438 | 0.0000000 |
| 8105 | GOCC MANNOSYLTRANSFERASE COMPLEX | 5 | 0.0026413 | 0.7763380 | 0.0073386 |
| 4563 | GOBP POSITIVE REGULATION OF NUCLEAR TRANSCRIBED MRNA CATABOLIC PROCESS DEADENYLATION DEPENDENT DECAY | 12 | 0.0000032 | 0.7760248 | 0.0000168 |
| 9970 | GOMF PROMOTER ENHANCER LOOP ANCHORING ACTIVITY | 6 | 0.0010602 | 0.7717184 | 0.0032924 |
| 8597 | GOCC U7 SNRNP | 7 | 0.0004097 | 0.7711505 | 0.0014179 |
| 8399 | GOCC PYRUVATE DEHYDROGENASE COMPLEX | 6 | 0.0010747 | 0.7708171 | 0.0033283 |
| 10048 | GOMF PYRIMIDINE NUCLEOTIDE BINDING | 5 | 0.0028463 | 0.7704504 | 0.0078317 |
head(subset(mres9$enrichment_result,p.adjustANOVA<0.05 & s.dist<0),50) |>
kbl(caption = "Top downregulated GO terms based on effect size (FDR<0.05) for H3K9ac") |>
kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 9376 | GOMF INHIBITORY MHC CLASS I RECEPTOR ACTIVITY | 8 | 0.0000039 | -0.9423883 | 0.0000200 |
| 8036 | GOCC IGG IMMUNOGLOBULIN COMPLEX | 13 | 0.0000003 | -0.8236985 | 0.0000016 |
| 1556 | GOBP FLAVONOID GLUCURONIDATION | 5 | 0.0032500 | -0.7599121 | 0.0088030 |
| 1531 | GOBP FC RECEPTOR MEDIATED INHIBITORY SIGNALING PATHWAY | 5 | 0.0039843 | -0.7434795 | 0.0105296 |
| 3837 | GOBP PEPTIDOGLYCAN METABOLIC PROCESS | 6 | 0.0024404 | -0.7143509 | 0.0068467 |
| 82 | GOBP ADENYLATE CYCLASE INHIBITING G PROTEIN COUPLED GLUTAMATE RECEPTOR SIGNALING PATHWAY | 9 | 0.0002185 | -0.7114465 | 0.0008087 |
| 9267 | GOMF G PROTEIN COUPLED GLUTAMATE RECEPTOR ACTIVITY | 9 | 0.0002185 | -0.7114465 | 0.0008087 |
| 8034 | GOCC IGD IMMUNOGLOBULIN COMPLEX | 6 | 0.0027810 | -0.7049999 | 0.0076783 |
| 8035 | GOCC IGE IMMUNOGLOBULIN COMPLEX | 6 | 0.0027886 | -0.7048028 | 0.0076953 |
| 9093 | GOMF DOPAMINE BINDING | 7 | 0.0013083 | -0.7014130 | 0.0039751 |
| 8727 | GOMF ALKANE 1 MONOOXYGENASE ACTIVITY | 5 | 0.0073758 | -0.6918361 | 0.0178385 |
| 8040 | GOCC IMMUNOGLOBULIN COMPLEX CIRCULATING | 12 | 0.0000540 | -0.6729651 | 0.0002256 |
| 7583 | GOBP XENOBIOTIC GLUCURONIDATION | 7 | 0.0021922 | -0.6684009 | 0.0062479 |
| 9737 | GOMF N FORMYL PEPTIDE RECEPTOR ACTIVITY | 5 | 0.0100157 | -0.6649599 | 0.0232422 |
| 1555 | GOBP FLAVONE METABOLIC PROCESS | 6 | 0.0055756 | -0.6533331 | 0.0140869 |
| 10150 | GOMF SIALIC ACID BINDING | 22 | 0.0000001 | -0.6512329 | 0.0000008 |
| 10285 | GOMF TRACE AMINE RECEPTOR ACTIVITY | 7 | 0.0030493 | -0.6465473 | 0.0083287 |
| 9572 | GOMF MHC CLASS I RECEPTOR ACTIVITY | 12 | 0.0001462 | -0.6329548 | 0.0005607 |
| 9544 | GOMF MELANOCORTIN RECEPTOR ACTIVITY | 6 | 0.0073005 | -0.6323723 | 0.0176973 |
| 8916 | GOMF CCR1 CHEMOKINE RECEPTOR BINDING | 5 | 0.0176196 | -0.6128635 | 0.0372271 |
| 8649 | GOMF 4 GALACTOSYL N ACETYLGLUCOSAMINIDE 3 ALPHA L FUCOSYLTRANSFERASE ACTIVITY | 6 | 0.0101310 | -0.6060939 | 0.0234472 |
| 8072 | GOCC KAINATE SELECTIVE GLUTAMATE RECEPTOR COMPLEX | 5 | 0.0190455 | -0.6054077 | 0.0397168 |
| 9437 | GOMF KAINATE SELECTIVE GLUTAMATE RECEPTOR ACTIVITY | 5 | 0.0190455 | -0.6054077 | 0.0397168 |
| 8032 | GOCC IGA IMMUNOGLOBULIN COMPLEX | 10 | 0.0009712 | -0.6022985 | 0.0030414 |
| 9216 | GOMF GLUTAMATE GATED CALCIUM ION CHANNEL ACTIVITY | 8 | 0.0048965 | -0.5743891 | 0.0126227 |
| 8074 | GOCC KERATIN FILAMENT | 92 | 0.0000000 | -0.5720903 | 0.0000000 |
| 10219 | GOMF STRUCTURAL CONSTITUENT OF SKIN EPIDERMIS | 37 | 0.0000000 | -0.5664672 | 0.0000000 |
| 3945 | GOBP PH REDUCTION | 6 | 0.0173132 | -0.5609953 | 0.0366838 |
| 3898 | GOBP PHENYLPROPANOID METABOLIC PROCESS | 7 | 0.0103938 | -0.5592005 | 0.0239651 |
| 9218 | GOMF GLUTAMATE RECEPTOR ACTIVITY | 26 | 0.0000009 | -0.5570584 | 0.0000049 |
| 2660 | GOBP MONOTERPENOID METABOLIC PROCESS | 6 | 0.0219683 | -0.5400007 | 0.0448525 |
| 9791 | GOMF OXIDOREDUCTASE ACTIVITY ACTING ON PAIRED DONORS WITH INCORPORATION OR REDUCTION OF MOLECULAR OXYGEN REDUCED IRON SULFUR PROTEIN AS ONE DONOR AND INCORPORATION OF ONE ATOM OF OXYGEN | 9 | 0.0051732 | -0.5381334 | 0.0132422 |
| 5793 | GOBP REGULATION OF INTESTINAL LIPID ABSORPTION | 10 | 0.0036150 | -0.5313165 | 0.0096514 |
| 5726 | GOBP REGULATION OF GUANYLATE CYCLASE ACTIVITY | 7 | 0.0150480 | -0.5305653 | 0.0325993 |
| 8873 | GOMF CALCIUM SENSITIVE GUANYLATE CYCLASE ACTIVATOR ACTIVITY | 6 | 0.0247486 | -0.5292499 | 0.0497213 |
| 3569 | GOBP NEUROTRANSMITTER LOADING INTO SYNAPTIC VESICLE | 7 | 0.0157758 | -0.5268209 | 0.0339153 |
| 5334 | GOBP REGULATION OF ATRIAL CARDIAC MUSCLE CELL MEMBRANE REPOLARIZATION | 8 | 0.0104272 | -0.5228613 | 0.0240314 |
| 9254 | GOMF GUANYLATE CYCLASE REGULATOR ACTIVITY | 9 | 0.0078300 | -0.5118198 | 0.0187413 |
| 8039 | GOCC IMMUNOGLOBULIN COMPLEX | 147 | 0.0000000 | -0.5109430 | 0.0000000 |
| 2100 | GOBP KERATINIZATION | 83 | 0.0000000 | -0.5071794 | 0.0000000 |
| 9691 | GOMF NMDA GLUTAMATE RECEPTOR ACTIVITY | 7 | 0.0204331 | -0.5059233 | 0.0421138 |
| 9839 | GOMF PEPTIDOGLYCAN MURALYTIC ACTIVITY | 14 | 0.0010678 | -0.5049304 | 0.0033111 |
| 8783 | GOMF AROMATASE ACTIVITY | 26 | 0.0000083 | -0.5048413 | 0.0000398 |
| 4189 | GOBP POSITIVE REGULATION OF COMPLEMENT ACTIVATION | 7 | 0.0214896 | -0.5017708 | 0.0439958 |
| 3926 | GOBP PHOSPHOLIPASE C ACTIVATING G PROTEIN COUPLED ACETYLCHOLINE RECEPTOR SIGNALING PATHWAY | 8 | 0.0145273 | -0.4989015 | 0.0316686 |
| 9217 | GOMF GLUTAMATE GATED RECEPTOR ACTIVITY | 16 | 0.0005878 | -0.4961695 | 0.0019538 |
| 81 | GOBP ADENYLATE CYCLASE INHIBITING G PROTEIN COUPLED ACETYLCHOLINE RECEPTOR SIGNALING PATHWAY | 8 | 0.0160029 | -0.4917360 | 0.0343414 |
| 9284 | GOMF HEMOGLOBIN BINDING | 9 | 0.0112454 | -0.4878910 | 0.0256116 |
| 9374 | GOMF IMMUNOGLOBULIN RECEPTOR BINDING | 19 | 0.0002582 | -0.4840449 | 0.0009384 |
| 1396 | GOBP EPOXYGENASE P450 PATHWAY | 18 | 0.0003774 | -0.4838842 | 0.0013180 |
par(mar=c(5,27,3,3))
top <- mres9$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)## [1] 5193up <- 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="H3K9ac",xlab="ES")
grid()
pdf("k9a_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="H3K9ac",xlab="ES")
grid()
dev.off()## png
## 2par(mar=c(5.1, 4.1, 4.1, 2.1))
if (! file.exists("k9a_mitch.html") ) {
mitch_report(mres9,"k9a_mitch.html")
}H3K36ac
m36 <- mitch_import(x=dk36,DEtype="deseq2",geneTable=gt2)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 92792## Note: no. genes in output = 59007## Note: estimated proportion of input genes in output = 0.636head(m36)## x
## 5_8S_rRNA -3.64204468
## 5S_rRNA -0.09835988
## 7SK 0.56711205
## A1BG -2.26646709
## A1BG-AS1 -2.44703622
## A1CF 1.69473867mres36 <- mitch_calc(x=m36,genesets=pw,minsetsize=5,cores=8,priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(subset(mres36$enrichment_result,p.adjustANOVA<0.05 & s.dist>0),50) |>
kbl(caption = "Top upregulated GO terms based on effect size (FDR<0.05) for H3K36ac") |>
kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 8556 | GOCC TRANSCRIPTION FACTOR AP 1 COMPLEX | 5 | 0.0094192 | 0.6704247 | 0.0487910 |
| 1099 | GOBP DETECTION OF MUSCLE STRETCH | 7 | 0.0072701 | 0.5857724 | 0.0398632 |
| 10353 | GOMF TYPE I INTERFERON RECEPTOR BINDING | 17 | 0.0000393 | 0.5757965 | 0.0005356 |
| 9739 | GOMF ODORANT BINDING | 119 | 0.0000000 | 0.4304729 | 0.0000000 |
| 1805 | GOBP HEART TRABECULA FORMATION | 15 | 0.0047780 | 0.4206717 | 0.0289711 |
| 9741 | GOMF OLFACTORY RECEPTOR ACTIVITY | 411 | 0.0000000 | 0.3909698 | 0.0000000 |
| 6513 | GOBP RELAXATION OF CARDIAC MUSCLE | 17 | 0.0056782 | 0.3873418 | 0.0330828 |
| 559 | GOBP CELLULAR GLUCURONIDATION | 21 | 0.0031980 | 0.3714780 | 0.0210024 |
| 6949 | GOBP SERINE PHOSPHORYLATION OF STAT PROTEIN | 23 | 0.0022865 | 0.3672713 | 0.0162220 |
| 6930 | GOBP SENSORY PERCEPTION OF SMELL | 438 | 0.0000000 | 0.3462336 | 0.0000000 |
| 6706 | GOBP RESPONSE TO MUSCLE STRETCH | 25 | 0.0051636 | 0.3229935 | 0.0306845 |
| 5096 | GOBP PROTEIN REFOLDING | 26 | 0.0058039 | 0.3124220 | 0.0336085 |
| 8190 | GOCC M BAND | 28 | 0.0052040 | 0.3049330 | 0.0308718 |
| 6926 | GOBP SENSORY PERCEPTION OF CHEMICAL STIMULUS | 516 | 0.0000000 | 0.2775635 | 0.0000000 |
| 1102 | GOBP DETECTION OF STIMULUS INVOLVED IN SENSORY PERCEPTION | 533 | 0.0000000 | 0.2702239 | 0.0000000 |
| 7664 | GOCC A BAND | 40 | 0.0030974 | 0.2701053 | 0.0205745 |
| 8585 | GOCC T CELL RECEPTOR COMPLEX | 143 | 0.0000121 | 0.2115168 | 0.0001919 |
| 1101 | GOBP DETECTION OF STIMULUS | 654 | 0.0000000 | 0.2014026 | 0.0000000 |
| 8068 | GOCC I BAND | 143 | 0.0048174 | 0.1362660 | 0.0291419 |
| 9271 | GOMF G PROTEIN COUPLED RECEPTOR ACTIVITY | 847 | 0.0000001 | 0.1058868 | 0.0000029 |
| 7803 | GOCC CONTRACTILE MUSCLE FIBER | 248 | 0.0065231 | 0.0999387 | 0.0368125 |
| 6924 | GOBP SENSORY PERCEPTION | 965 | 0.0001135 | 0.0723280 | 0.0013496 |
head(subset(mres9$enrichment_result,p.adjustANOVA<0.05 & s.dist<0),50) |>
kbl(caption = "Top downregulated GO terms based on effect size (FDR<0.05) for H3K36ac") |>
kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 9376 | GOMF INHIBITORY MHC CLASS I RECEPTOR ACTIVITY | 8 | 0.0000039 | -0.9423883 | 0.0000200 |
| 8036 | GOCC IGG IMMUNOGLOBULIN COMPLEX | 13 | 0.0000003 | -0.8236985 | 0.0000016 |
| 1556 | GOBP FLAVONOID GLUCURONIDATION | 5 | 0.0032500 | -0.7599121 | 0.0088030 |
| 1531 | GOBP FC RECEPTOR MEDIATED INHIBITORY SIGNALING PATHWAY | 5 | 0.0039843 | -0.7434795 | 0.0105296 |
| 3837 | GOBP PEPTIDOGLYCAN METABOLIC PROCESS | 6 | 0.0024404 | -0.7143509 | 0.0068467 |
| 82 | GOBP ADENYLATE CYCLASE INHIBITING G PROTEIN COUPLED GLUTAMATE RECEPTOR SIGNALING PATHWAY | 9 | 0.0002185 | -0.7114465 | 0.0008087 |
| 9267 | GOMF G PROTEIN COUPLED GLUTAMATE RECEPTOR ACTIVITY | 9 | 0.0002185 | -0.7114465 | 0.0008087 |
| 8034 | GOCC IGD IMMUNOGLOBULIN COMPLEX | 6 | 0.0027810 | -0.7049999 | 0.0076783 |
| 8035 | GOCC IGE IMMUNOGLOBULIN COMPLEX | 6 | 0.0027886 | -0.7048028 | 0.0076953 |
| 9093 | GOMF DOPAMINE BINDING | 7 | 0.0013083 | -0.7014130 | 0.0039751 |
| 8727 | GOMF ALKANE 1 MONOOXYGENASE ACTIVITY | 5 | 0.0073758 | -0.6918361 | 0.0178385 |
| 8040 | GOCC IMMUNOGLOBULIN COMPLEX CIRCULATING | 12 | 0.0000540 | -0.6729651 | 0.0002256 |
| 7583 | GOBP XENOBIOTIC GLUCURONIDATION | 7 | 0.0021922 | -0.6684009 | 0.0062479 |
| 9737 | GOMF N FORMYL PEPTIDE RECEPTOR ACTIVITY | 5 | 0.0100157 | -0.6649599 | 0.0232422 |
| 1555 | GOBP FLAVONE METABOLIC PROCESS | 6 | 0.0055756 | -0.6533331 | 0.0140869 |
| 10150 | GOMF SIALIC ACID BINDING | 22 | 0.0000001 | -0.6512329 | 0.0000008 |
| 10285 | GOMF TRACE AMINE RECEPTOR ACTIVITY | 7 | 0.0030493 | -0.6465473 | 0.0083287 |
| 9572 | GOMF MHC CLASS I RECEPTOR ACTIVITY | 12 | 0.0001462 | -0.6329548 | 0.0005607 |
| 9544 | GOMF MELANOCORTIN RECEPTOR ACTIVITY | 6 | 0.0073005 | -0.6323723 | 0.0176973 |
| 8916 | GOMF CCR1 CHEMOKINE RECEPTOR BINDING | 5 | 0.0176196 | -0.6128635 | 0.0372271 |
| 8649 | GOMF 4 GALACTOSYL N ACETYLGLUCOSAMINIDE 3 ALPHA L FUCOSYLTRANSFERASE ACTIVITY | 6 | 0.0101310 | -0.6060939 | 0.0234472 |
| 8072 | GOCC KAINATE SELECTIVE GLUTAMATE RECEPTOR COMPLEX | 5 | 0.0190455 | -0.6054077 | 0.0397168 |
| 9437 | GOMF KAINATE SELECTIVE GLUTAMATE RECEPTOR ACTIVITY | 5 | 0.0190455 | -0.6054077 | 0.0397168 |
| 8032 | GOCC IGA IMMUNOGLOBULIN COMPLEX | 10 | 0.0009712 | -0.6022985 | 0.0030414 |
| 9216 | GOMF GLUTAMATE GATED CALCIUM ION CHANNEL ACTIVITY | 8 | 0.0048965 | -0.5743891 | 0.0126227 |
| 8074 | GOCC KERATIN FILAMENT | 92 | 0.0000000 | -0.5720903 | 0.0000000 |
| 10219 | GOMF STRUCTURAL CONSTITUENT OF SKIN EPIDERMIS | 37 | 0.0000000 | -0.5664672 | 0.0000000 |
| 3945 | GOBP PH REDUCTION | 6 | 0.0173132 | -0.5609953 | 0.0366838 |
| 3898 | GOBP PHENYLPROPANOID METABOLIC PROCESS | 7 | 0.0103938 | -0.5592005 | 0.0239651 |
| 9218 | GOMF GLUTAMATE RECEPTOR ACTIVITY | 26 | 0.0000009 | -0.5570584 | 0.0000049 |
| 2660 | GOBP MONOTERPENOID METABOLIC PROCESS | 6 | 0.0219683 | -0.5400007 | 0.0448525 |
| 9791 | GOMF OXIDOREDUCTASE ACTIVITY ACTING ON PAIRED DONORS WITH INCORPORATION OR REDUCTION OF MOLECULAR OXYGEN REDUCED IRON SULFUR PROTEIN AS ONE DONOR AND INCORPORATION OF ONE ATOM OF OXYGEN | 9 | 0.0051732 | -0.5381334 | 0.0132422 |
| 5793 | GOBP REGULATION OF INTESTINAL LIPID ABSORPTION | 10 | 0.0036150 | -0.5313165 | 0.0096514 |
| 5726 | GOBP REGULATION OF GUANYLATE CYCLASE ACTIVITY | 7 | 0.0150480 | -0.5305653 | 0.0325993 |
| 8873 | GOMF CALCIUM SENSITIVE GUANYLATE CYCLASE ACTIVATOR ACTIVITY | 6 | 0.0247486 | -0.5292499 | 0.0497213 |
| 3569 | GOBP NEUROTRANSMITTER LOADING INTO SYNAPTIC VESICLE | 7 | 0.0157758 | -0.5268209 | 0.0339153 |
| 5334 | GOBP REGULATION OF ATRIAL CARDIAC MUSCLE CELL MEMBRANE REPOLARIZATION | 8 | 0.0104272 | -0.5228613 | 0.0240314 |
| 9254 | GOMF GUANYLATE CYCLASE REGULATOR ACTIVITY | 9 | 0.0078300 | -0.5118198 | 0.0187413 |
| 8039 | GOCC IMMUNOGLOBULIN COMPLEX | 147 | 0.0000000 | -0.5109430 | 0.0000000 |
| 2100 | GOBP KERATINIZATION | 83 | 0.0000000 | -0.5071794 | 0.0000000 |
| 9691 | GOMF NMDA GLUTAMATE RECEPTOR ACTIVITY | 7 | 0.0204331 | -0.5059233 | 0.0421138 |
| 9839 | GOMF PEPTIDOGLYCAN MURALYTIC ACTIVITY | 14 | 0.0010678 | -0.5049304 | 0.0033111 |
| 8783 | GOMF AROMATASE ACTIVITY | 26 | 0.0000083 | -0.5048413 | 0.0000398 |
| 4189 | GOBP POSITIVE REGULATION OF COMPLEMENT ACTIVATION | 7 | 0.0214896 | -0.5017708 | 0.0439958 |
| 3926 | GOBP PHOSPHOLIPASE C ACTIVATING G PROTEIN COUPLED ACETYLCHOLINE RECEPTOR SIGNALING PATHWAY | 8 | 0.0145273 | -0.4989015 | 0.0316686 |
| 9217 | GOMF GLUTAMATE GATED RECEPTOR ACTIVITY | 16 | 0.0005878 | -0.4961695 | 0.0019538 |
| 81 | GOBP ADENYLATE CYCLASE INHIBITING G PROTEIN COUPLED ACETYLCHOLINE RECEPTOR SIGNALING PATHWAY | 8 | 0.0160029 | -0.4917360 | 0.0343414 |
| 9284 | GOMF HEMOGLOBIN BINDING | 9 | 0.0112454 | -0.4878910 | 0.0256116 |
| 9374 | GOMF IMMUNOGLOBULIN RECEPTOR BINDING | 19 | 0.0002582 | -0.4840449 | 0.0009384 |
| 1396 | GOBP EPOXYGENASE P450 PATHWAY | 18 | 0.0003774 | -0.4838842 | 0.0013180 |
par(mar=c(5,27,3,3))
top <- mres36$enrichment_result
top <- subset(top,p.adjustANOVA<0.05)
nrow(top)## [1] 2041up <- 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))
names(vec) <- stringr::str_trunc(names(vec), 70)
vec <- vec[order(vec)]
barplot(abs(vec),col=sign(-vec)+3,horiz=TRUE,las=1,cex.names=0.65,main="H3K36ac",xlab="ES")
grid()
pdf("k36a_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="H3K36ac",xlab="ES")
grid()
dev.off()## png
## 2par(mar=c(5.1, 4.1, 4.1, 2.1))
if (! file.exists("k36a_mitch.html") ) {
mitch_report(mres36,"k36a_mitch.html")
}SLC
rpm9_SLC2A4 <- rpm9[grep("SLC2A4",rownames(rpm9)),]
rpm9_SLC2A4## pos1 pre1 pos2
## ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 26.70826 19.95114 26.00418
## ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 37.49235 17.23053 33.11279
## pre2 pos3 pre3
## ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 20.57199 27.83474 20.59174
## ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 16.66161 39.55999 17.82569
## pos4 pre4 pos5
## ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 19.57435 32.65682 14.91572
## ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 36.32246 19.81503 18.37829
## pre5 pos6 pre6
## ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 24.01717 25.25616 29.76552
## ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 19.02484 33.77164 36.26257rpm9_SLC2A4 <- colSums(rpm9_SLC2A4)
rpm9_SLC2A4 <- matrix(rpm9_SLC2A4,nrow=2)
colnames(rpm9_SLC2A4) <- 1:6
rownames(rpm9_SLC2A4) <- c("Post","Pre")
rpm9_SLC2A4 <- rpm9_SLC2A4[c(2,1),]
barplot(rpm9_SLC2A4,beside=TRUE,col=c("#0047AB", "#D2042D"),
main="SLC2A4 H3K9ac",ylab="Counts per million",
xlab="Participant",ylim=c(0,73))
legend(11.5,73, legend=c("Pre", "Post"), fill=c("#0047AB", "#D2042D"), cex=1)
pdf("k9a_SLC2A4.pdf",width=5,height=5)
barplot(rpm9_SLC2A4,beside=TRUE,col=c("#0047AB", "#D2042D"),
main="SLC2A4 H3K9ac",ylab="Counts per million",
xlab="Participant",ylim=c(0,73))
legend(11.5,73, legend=c("Pre", "Post"), fill=c("#0047AB", "#D2042D"), cex=1)
dev.off()## png
## 2rpm36_SLC2A4 <- rpm36[grep("SLC2A4",rownames(rpm36)),]
rpm36_SLC2A4## pos1 pre1
## ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 17.01140 17.783909
## ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 32.22484 7.948869
## pos2 pre2
## ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 9.363935 19.377948
## ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 22.719057 9.527491
## pos3 pre3 pos4
## ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 15.24830 19.12591 11.42846
## ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 50.48225 19.25427 39.13004
## pre4 pos5 pre5
## ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 20.06859 21.25035 24.76518
## ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 17.42318 46.66743 14.45623
## pos6 pre6
## ENSG00000125520.14 SLC2A4RG chr20 63737776 63744528 25.28059 24.06254
## ENSG00000181856.15 SLC2A4 chr17 7279718 7284427 34.14333 20.83766rpm36_SLC2A4 <- colSums(rpm36_SLC2A4)
rpm36_SLC2A4 <- matrix(rpm36_SLC2A4,nrow=2)
colnames(rpm36_SLC2A4) <- 1:6
rownames(rpm36_SLC2A4) <- c("Post","Pre")
rpm36_SLC2A4 <- rpm36_SLC2A4[c(2,1),]
barplot(rpm36_SLC2A4,beside=TRUE,col=c("#0047AB", "#D2042D"),
main="SLC2A4 H3K36ac",ylab="counts per million",
xlab="participant",ylim=c(0,70))
legend(3.5,70, legend=c("Pre", "Post"), fill=c("#0047AB","#D2042D"), cex=1)
pdf("k36a_SLC2A4.pdf",width=5,height=5)
barplot(rpm36_SLC2A4,beside=TRUE,col=c("#0047AB", "#D2042D"),
main="SLC2A4 H3K36ac",ylab="counts per million",
xlab="participant",ylim=c(0,70))
legend(3.5,70, legend=c("Pre", "Post"), fill=c("#0047AB","#D2042D"), cex=1)
dev.off()## png
## 2Session
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/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 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: Etc/UTC
## 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.4
## [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.1 GenomeInfoDb_1.40.1
## [15] IRanges_2.38.1 S4Vectors_0.42.1
## [17] BiocGenerics_0.50.0 reshape2_1.4.4
## [19] RhpcBLASctl_0.23-42
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.2 dplyr_1.1.4
## [4] bitops_1.0-8 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.8
## [19] knitr_1.47 S4Arrays_1.4.1 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-1 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.8
## [82] xfun_0.45 pkgconfig_2.0.3