ECHS1 KO RNA-seq
Mark Ziemann
2022-01-19
Source: https://github.com/markziemann/echs1_rnaseq
Introduction
The samples are labelled as CON (control) and KO for knockout, ut are the untreated samples and dNs are the treated samples (mito biogenesis stimulator).
suppressPackageStartupMessages({
library("zoo")
library("tidyverse")
library("reshape2")
library("DESeq2")
library("gplots")
library("fgsea")
library("MASS")
library("mitch")
library("eulerr")
library("limma")
library("topconfects")
library("kableExtra")
})Import read counts
Importing RNA-seq data
tmp <- read.table("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
xx$geneid = NULL
xx <- round(xx)
xx <- xx[,which(colnames(xx)!="test")]
xx[1:6,1:6]## 1_CON_ut_A_RNA_HB_HLF57DRXY_CTGTTGGTCC-AACGGTCTAT_L001
## ENSG00000000003.15 TSPAN6 704
## ENSG00000000005.6 TNMD 0
## ENSG00000000419.13 DPM1 624
## ENSG00000000457.14 SCYL3 97
## ENSG00000000460.17 C1orf112 110
## ENSG00000000938.13 FGR 0
## 1_CON_ut_A_RNA_HB_HLF57DRXY_CTGTTGGTCC-AACGGTCTAT_L002
## ENSG00000000003.15 TSPAN6 676
## ENSG00000000005.6 TNMD 0
## ENSG00000000419.13 DPM1 643
## ENSG00000000457.14 SCYL3 99
## ENSG00000000460.17 C1orf112 122
## ENSG00000000938.13 FGR 0
## 10_KO_ut_B_RNA_HB_HLF57DRXY_CACGGAACAA-GTGCTAGGTT_L001
## ENSG00000000003.15 TSPAN6 752
## ENSG00000000005.6 TNMD 0
## ENSG00000000419.13 DPM1 360
## ENSG00000000457.14 SCYL3 102
## ENSG00000000460.17 C1orf112 110
## ENSG00000000938.13 FGR 0
## 10_KO_ut_B_RNA_HB_HLF57DRXY_CACGGAACAA-GTGCTAGGTT_L002
## ENSG00000000003.15 TSPAN6 803
## ENSG00000000005.6 TNMD 0
## ENSG00000000419.13 DPM1 334
## ENSG00000000457.14 SCYL3 75
## ENSG00000000460.17 C1orf112 100
## ENSG00000000938.13 FGR 0
## 11_KO_ut_C_RNA_HB_HLF57DRXY_TGGAGTACTT-TCCACACAGA_L001
## ENSG00000000003.15 TSPAN6 1502
## ENSG00000000005.6 TNMD 0
## ENSG00000000419.13 DPM1 703
## ENSG00000000457.14 SCYL3 203
## ENSG00000000460.17 C1orf112 175
## ENSG00000000938.13 FGR 0
## 11_KO_ut_C_RNA_HB_HLF57DRXY_TGGAGTACTT-TCCACACAGA_L002
## ENSG00000000003.15 TSPAN6 1509
## ENSG00000000005.6 TNMD 0
## ENSG00000000419.13 DPM1 659
## ENSG00000000457.14 SCYL3 184
## ENSG00000000460.17 C1orf112 191
## ENSG00000000938.13 FGR 0dim(xx)## [1] 60651 32Fix the sample names.
They are duplicated for lane 1 and 2, which I will aggregate.
colnames(xx) <- sapply(strsplit(colnames(xx),"_RNA_"),"[[",1)
colnames(xx) <- sapply(strsplit(sub("_","@",colnames(xx)),"@"),"[[",2)
labels <- unique(sapply(strsplit(colnames(xx),"\\."),"[[",1))
l <- lapply(labels,function(x) { rowSums(xx[,grep(x,colnames(xx))]) } )
ll <- do.call(cbind,l)
colnames(ll) <- labels
ll <- as.data.frame(ll[,order(colnames(ll))])Make a sample sheet. This will be easy because the groups are specified in the name.
ss <- as.data.frame(colnames(ll))
ss$ko <- factor(grepl("KO", colnames(ll)))
ss$trt <- factor(grepl("dNs",colnames(ll)))
rownames(ss) <- ss[,1]
ss[,1]=NULL
# can also set some colours for later
ss$cols <- c(rep("gray",4),rep("lightblue",4),rep("pink",4),rep("orange",4))
ss %>% kbl(caption = "Sample sheet") %>% kable_paper("hover", full_width = F)| ko | trt | cols | |
|---|---|---|---|
| CON_dNs_A | FALSE | TRUE | gray |
| CON_dNs_B | FALSE | TRUE | gray |
| CON_dNs_C | FALSE | TRUE | gray |
| CON_dNs_D | FALSE | TRUE | gray |
| CON_ut_A | FALSE | FALSE | lightblue |
| CON_ut_B | FALSE | FALSE | lightblue |
| CON_ut_C | FALSE | FALSE | lightblue |
| CON_ut_D | FALSE | FALSE | lightblue |
| KO_dNs_A | TRUE | TRUE | pink |
| KO_dNs_B | TRUE | TRUE | pink |
| KO_dNs_C | TRUE | TRUE | pink |
| KO_dNs_D | TRUE | TRUE | pink |
| KO_ut_A | TRUE | FALSE | orange |
| KO_ut_B | TRUE | FALSE | orange |
| KO_ut_C | TRUE | FALSE | orange |
| KO_ut_D | TRUE | FALSE | orange |
QC analysis
Here I’ll look at a few different quality control measures.
par(mar=c(5,8,3,1))
barplot(colSums(ll),horiz=TRUE,las=1,xlab="num reads",col=ss$cols)
sums <- colSums(ll)
sums <- sums[order(sums)]
barplot(sums,horiz=TRUE,las=1,xlab="num reads")
abline(v=15000000,col="red")
MDS plot for all samples
Multidimensional scaling plot to show the variation between all samples, very similar to PCA.
Firstly with the data before aggregating technical replicates.
mds <- cmdscale(dist(t(xx)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",bty="n",pch=19, cex=4 ,col="gray")
text(mds, labels=rownames(mds) )
Now for the data after aggregation.
It looks like the drug has a bigger effect than the gene KO.
mds <- cmdscale(dist(t(ll)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",bty="n",pch=19, cex=4 ,col=ss$cols)
text(mds, labels=rownames(mds) )
Correlation heatmap
heatmap.2(cor(ll),trace="n",main="Pearson correlation heatmap",margin=c(8,8))
MDS plot for UT samples
ll2 <- ll[,grep("ut",colnames(ll))]
mds <- cmdscale(dist(t(ll2)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",bty="n",pch=19, cex=4 ,col=ss$cols)
text(mds, labels=rownames(mds) )
MDS plot for dNs samples
ll3 <- ll[,grep("dNs",colnames(ll))]
mds <- cmdscale(dist(t(ll3)))
plot(mds, xlab="Coordinate 1", ylab="Coordinate 2",
type = "p",bty="n",pch=19, cex=4 ,col=ss$cols)
text(mds, labels=rownames(mds) )
Correlation heatmaps for UT and dNs separately
heatmap.2(cor(ll2),trace="n",main="Pearson correlation heatmap",margin=c(8,8))
heatmap.2(cor(ll3),trace="n",main="Pearson correlation heatmap",margin=c(8,8))
Set up the different datasets for differential expression analysis
Don’t forget to remove poorly detected genes from the matrix with a threshold of 10 reads per sample on average.
There are 4 contrasts to set up.
Effect of the KO in untreated cells
Effect of the KO in treated cells
Effect of the drug in WT cells
4, Effect of the drug in KO cells
ss1 <- ss[which(ss$trt=="FALSE"),]
xx1 <- ll[which(colnames(ll) %in% rownames(ss1))]
xx1 <- xx1[which(rowSums(xx1)>=10),]
ss2 <- ss[which(ss$trt=="TRUE"),]
xx2 <- ll[which(colnames(ll) %in% rownames(ss2))]
xx2 <- xx2[which(rowSums(xx2)>=10),]
ss3 <- ss[which(ss$ko=="FALSE"),]
xx3 <- ll[which(colnames(ll) %in% rownames(ss3))]
xx3 <- xx3[which(rowSums(xx3)>=10),]
ss4 <- ss[which(ss$ko=="TRUE"),]
xx4 <- ll[which(colnames(ll) %in% rownames(ss4))]
xx4 <- xx4[which(rowSums(xx4)>=10),]Differential expression with DESeq2
dds <- DESeqDataSetFromMatrix(countData = xx1 , colData = ss1, design = ~ ko )## 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 <- as.data.frame(zz[order(zz$pvalue),])
dge[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 1: Effect of the KO in untreated cells") %>% kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| ENSG00000049192.15 ADAMTS6 | 1232.5643 | -3.320315 | 0.0735209 | -45.16150 | 0 | 0 |
| ENSG00000090339.9 ICAM1 | 8426.0162 | 2.137068 | 0.0421343 | 50.72040 | 0 | 0 |
| ENSG00000099875.16 MKNK2 | 5196.7000 | 1.646797 | 0.0435241 | 37.83644 | 0 | 0 |
| ENSG00000105825.14 TFPI2 | 10749.0439 | -2.758378 | 0.0645077 | -42.76045 | 0 | 0 |
| ENSG00000107551.21 RASSF4 | 1224.3947 | 3.078164 | 0.0812402 | 37.88965 | 0 | 0 |
| ENSG00000116016.14 EPAS1 | 16685.0850 | 1.887592 | 0.0496555 | 38.01373 | 0 | 0 |
| ENSG00000125538.12 IL1B | 9573.0006 | -1.925913 | 0.0457459 | -42.10027 | 0 | 0 |
| ENSG00000138080.14 EMILIN1 | 6299.0508 | 2.088516 | 0.0483164 | 43.22583 | 0 | 0 |
| ENSG00000146197.9 SCUBE3 | 7326.2924 | -2.787733 | 0.0674659 | -41.32061 | 0 | 0 |
| ENSG00000154146.13 NRGN | 1528.8050 | 2.283661 | 0.0600431 | 38.03371 | 0 | 0 |
| ENSG00000164949.8 GEM | 1954.1763 | -2.372324 | 0.0565928 | -41.91920 | 0 | 0 |
| ENSG00000165507.9 DEPP1 | 3368.7544 | 2.850455 | 0.0522246 | 54.58071 | 0 | 0 |
| ENSG00000166250.12 CLMP | 20860.4032 | -1.674644 | 0.0397766 | -42.10125 | 0 | 0 |
| ENSG00000166670.10 MMP10 | 3155.5769 | -6.120693 | 0.1021577 | -59.91416 | 0 | 0 |
| ENSG00000166920.13 C15orf48 | 2293.5579 | 2.687532 | 0.0645990 | 41.60330 | 0 | 0 |
| ENSG00000167772.12 ANGPTL4 | 14717.2697 | 3.591520 | 0.0856872 | 41.91430 | 0 | 0 |
| ENSG00000168528.12 SERINC2 | 4459.4824 | 1.937358 | 0.0495824 | 39.07346 | 0 | 0 |
| ENSG00000170961.7 HAS2 | 3615.4379 | -4.322465 | 0.0852007 | -50.73275 | 0 | 0 |
| ENSG00000172572.7 PDE3A | 1574.2471 | -3.854290 | 0.0926037 | -41.62136 | 0 | 0 |
| ENSG00000172927.8 MYEOV | 928.0298 | 4.057927 | 0.0905625 | 44.80804 | 0 | 0 |
dge1 <- dge
d1up <- rownames(subset(dge1,padj <= 0.05 & log2FoldChange > 0))
d1dn <- rownames(subset(dge1,padj <= 0.05 & log2FoldChange < 0))
write.table(dge1,file="dge1.tsv",quote=FALSE,sep="\t")dds <- DESeqDataSetFromMatrix(countData = xx2 , colData = ss2, design = ~ ko )## 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 <- as.data.frame(zz[order(zz$pvalue),])
dge[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 2: Effect of the KO in treated cells") %>% kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| ENSG00000081248.12 CACNA1S | 1209.4345 | 4.728860 | 0.1048201 | 45.11406 | 0 | 0 |
| ENSG00000142910.16 TINAGL1 | 4889.9014 | 3.128308 | 0.0629648 | 49.68348 | 0 | 0 |
| ENSG00000162804.14 SNED1 | 6126.7547 | -2.358713 | 0.0490839 | -48.05475 | 0 | 0 |
| ENSG00000166920.13 C15orf48 | 2299.7911 | 3.960911 | 0.0886313 | 44.68975 | 0 | 0 |
| ENSG00000167772.12 ANGPTL4 | 9936.5018 | 3.543783 | 0.0816955 | 43.37797 | 0 | 0 |
| ENSG00000113361.13 CDH6 | 816.5796 | -3.817982 | 0.1059909 | -36.02179 | 0 | 0 |
| ENSG00000101187.16 SLCO4A1 | 19101.5516 | 1.934484 | 0.0565000 | 34.23868 | 0 | 0 |
| ENSG00000143127.13 ITGA10 | 1156.7355 | -3.687434 | 0.1104110 | -33.39734 | 0 | 0 |
| ENSG00000151062.15 CACNA2D4 | 655.2444 | 4.551764 | 0.1378037 | 33.03079 | 0 | 0 |
| ENSG00000170961.7 HAS2 | 613.5916 | -3.363180 | 0.1045664 | -32.16312 | 0 | 0 |
| ENSG00000168528.12 SERINC2 | 3732.8640 | 1.667401 | 0.0525189 | 31.74856 | 0 | 0 |
| ENSG00000135245.10 HILPDA | 1510.8969 | 2.621093 | 0.0845048 | 31.01708 | 0 | 0 |
| ENSG00000172927.8 MYEOV | 472.9211 | 4.325653 | 0.1395499 | 30.99717 | 0 | 0 |
| ENSG00000115756.13 HPCAL1 | 4084.1250 | 1.608911 | 0.0524306 | 30.68650 | 0 | 0 |
| ENSG00000100985.7 MMP9 | 14795.4099 | -2.331039 | 0.0761929 | -30.59392 | 0 | 0 |
| ENSG00000163235.16 TGFA | 1660.7730 | 2.758967 | 0.0910679 | 30.29573 | 0 | 0 |
| ENSG00000125571.9 IL37 | 1230.5444 | 2.713549 | 0.0896194 | 30.27858 | 0 | 0 |
| ENSG00000118257.17 NRP2 | 2033.2565 | -2.000568 | 0.0681154 | -29.37026 | 0 | 0 |
| ENSG00000213658.12 LAT | 1708.1308 | -1.949041 | 0.0664827 | -29.31650 | 0 | 0 |
| ENSG00000101443.18 WFDC2 | 5495.9129 | 2.267593 | 0.0773716 | 29.30783 | 0 | 0 |
dge2 <- dge
d2up <- rownames(subset(dge2,padj <= 0.05 & log2FoldChange > 0))
d2dn <- rownames(subset(dge2,padj <= 0.05 & log2FoldChange < 0))
write.table(dge2,file="dge2.tsv",quote=FALSE,sep="\t")dds <- DESeqDataSetFromMatrix(countData = xx3 , colData = ss3, design = ~ trt )## 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 <- as.data.frame(zz[order(zz$pvalue),])
dge[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 3: Effect of the drug in WT cells") %>% kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| ENSG00000068912.14 ERLEC1 | 1851.848 | -2.762775 | 0.0670700 | -41.19240 | 0 | 0 |
| ENSG00000078747.16 ITCH | 2736.282 | -2.038664 | 0.0523301 | -38.95780 | 0 | 0 |
| ENSG00000100836.10 PABPN1 | 6539.865 | 2.072716 | 0.0487025 | 42.55874 | 0 | 0 |
| ENSG00000102755.12 FLT1 | 5426.297 | -3.705474 | 0.0880185 | -42.09882 | 0 | 0 |
| ENSG00000106366.9 SERPINE1 | 24625.790 | -2.091105 | 0.0507831 | -41.17717 | 0 | 0 |
| ENSG00000108821.14 COL1A1 | 9382.064 | 2.181043 | 0.0456926 | 47.73301 | 0 | 0 |
| ENSG00000110048.12 OSBP | 3405.995 | -2.355956 | 0.0614408 | -38.34517 | 0 | 0 |
| ENSG00000112473.18 SLC39A7 | 4521.830 | -2.281225 | 0.0578603 | -39.42640 | 0 | 0 |
| ENSG00000116285.13 ERRFI1 | 36944.117 | -2.684910 | 0.0567602 | -47.30267 | 0 | 0 |
| ENSG00000122884.12 P4HA1 | 4295.439 | -2.579782 | 0.0582809 | -44.26463 | 0 | 0 |
| ENSG00000131016.17 AKAP12 | 21452.521 | -3.317855 | 0.0663163 | -50.03076 | 0 | 0 |
| ENSG00000138166.6 DUSP5 | 4318.320 | -2.155073 | 0.0498547 | -43.22711 | 0 | 0 |
| ENSG00000142949.17 PTPRF | 7682.322 | -1.848612 | 0.0487213 | -37.94262 | 0 | 0 |
| ENSG00000150961.15 SEC24D | 2402.484 | -3.024974 | 0.0617986 | -48.94889 | 0 | 0 |
| ENSG00000150995.20 ITPR1 | 4106.314 | -3.021908 | 0.0665647 | -45.39805 | 0 | 0 |
| ENSG00000159167.12 STC1 | 13523.272 | -3.145427 | 0.0722957 | -43.50780 | 0 | 0 |
| ENSG00000162804.14 SNED1 | 6283.910 | 2.862292 | 0.0481808 | 59.40725 | 0 | 0 |
| ENSG00000164484.12 TMEM200A | 1197.746 | -3.470114 | 0.0881540 | -39.36421 | 0 | 0 |
| ENSG00000166250.12 CLMP | 18283.547 | -2.335014 | 0.0538393 | -43.37009 | 0 | 0 |
| ENSG00000166670.10 MMP10 | 3245.447 | -3.537259 | 0.0771167 | -45.86894 | 0 | 0 |
dge3 <- dge
d3up <- rownames(subset(dge3,padj <= 0.05 & log2FoldChange > 0))
d3dn <- rownames(subset(dge3,padj <= 0.05 & log2FoldChange < 0))
write.table(dge3,file="dge3.tsv",quote=FALSE,sep="\t")dds <- DESeqDataSetFromMatrix(countData = xx4 , colData = ss4, design = ~ trt )## 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 <- as.data.frame(zz[order(zz$pvalue),])
dge[1:20,1:6] %>% kbl(caption = "Top gene expression differences for contrast 3: Effect of the drug in KO cells") %>% kable_paper("hover", full_width = F)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| ENSG00000059804.16 SLC2A3 | 5623.218 | -1.769194 | 0.0468174 | -37.78920 | 0 | 0 |
| ENSG00000103888.17 CEMIP | 25230.452 | -2.601566 | 0.0673462 | -38.62975 | 0 | 0 |
| ENSG00000104419.17 NDRG1 | 38083.260 | -2.549397 | 0.0511771 | -49.81517 | 0 | 0 |
| ENSG00000134250.20 NOTCH2 | 6533.444 | -2.176937 | 0.0578058 | -37.65947 | 0 | 0 |
| ENSG00000197324.9 LRP10 | 9816.752 | -1.770453 | 0.0378724 | -46.74781 | 0 | 0 |
| ENSG00000198712.1 MT-CO2 | 290785.069 | 2.997493 | 0.0744604 | 40.25620 | 0 | 0 |
| ENSG00000198727.2 MT-CYB | 116486.352 | 3.700459 | 0.0765560 | 48.33662 | 0 | 0 |
| ENSG00000198763.3 MT-ND2 | 121199.785 | 3.262193 | 0.0678250 | 48.09722 | 0 | 0 |
| ENSG00000198840.2 MT-ND3 | 17981.096 | 2.878705 | 0.0705336 | 40.81322 | 0 | 0 |
| ENSG00000198959.12 TGM2 | 1712.100 | -2.516008 | 0.0594505 | -42.32104 | 0 | 0 |
| ENSG00000225630.1 MTND2P28 | 5023.069 | 3.268132 | 0.0682697 | 47.87088 | 0 | 0 |
| ENSG00000258017.2 AC011603.2 | 7033.925 | 2.918598 | 0.0767330 | 38.03578 | 0 | 0 |
| ENSG00000259207.9 ITGB3 | 6216.236 | -2.380112 | 0.0532424 | -44.70328 | 0 | 0 |
| ENSG00000198804.2 MT-CO1 | 506735.846 | 2.899638 | 0.0779490 | 37.19916 | 0 | 0 |
| ENSG00000134531.10 EMP1 | 8062.151 | -1.864264 | 0.0501448 | -37.17758 | 0 | 0 |
| ENSG00000198888.2 MT-ND1 | 70406.196 | 2.780834 | 0.0758103 | 36.68147 | 0 | 0 |
| ENSG00000142949.17 PTPRF | 7521.711 | -1.802689 | 0.0498077 | -36.19295 | 0 | 0 |
| ENSG00000198938.2 MT-CO3 | 211300.735 | 2.731621 | 0.0754794 | 36.19029 | 0 | 0 |
| ENSG00000178726.7 THBD | 1419.592 | -2.752535 | 0.0765192 | -35.97182 | 0 | 0 |
| ENSG00000131016.17 AKAP12 | 17614.587 | -2.494844 | 0.0707157 | -35.27991 | 0 | 0 |
dge4 <- dge
d4up <- rownames(subset(dge4,padj <= 0.05 & log2FoldChange > 0))
d4dn <- rownames(subset(dge4,padj <= 0.05 & log2FoldChange < 0))
write.table(dge4,file="dge4.tsv",quote=FALSE,sep="\t")Here let’s look at some plots.
MA plot shows the average level and fold change of all detected genes. Volcano plot shows the fold change and the significance, as measured by -log(p-value). Significant genes are shown as red points.
There are heatmaps of the top ranked genes by p-value. Above the gene expression values there is a bar in yellow/orange colours. Yellow shows relatively low QuickDASH score and orange shows high score.
maplot <- function(de,contrast_name) {
sig <-subset(de, padj < 0.05 )
up <-rownames(subset(de, padj < 0.05 & log2FoldChange > 0))
dn <-rownames(subset(de, padj < 0.05 & log2FoldChange < 0))
GENESUP <- length(up)
GENESDN <- length(dn)
DET=nrow(de)
SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down", DET, "detected")
ns <-subset(de, padj > 0.05 )
plot(log2(de$baseMean),de$log2FoldChange,
xlab="log2 basemean", ylab="log2 foldchange",
pch=19, cex=0.5, col="dark gray",
main=contrast_name, cex.main=1)
points(log2(sig$baseMean),sig$log2FoldChange,
pch=19, cex=0.5, col="red")
mtext(SUBHEADER,cex = 1)
}
make_volcano <- function(de,name) {
sig <- subset(de,padj<0.05)
N_SIG=nrow(sig)
N_UP=nrow(subset(sig,log2FoldChange>0))
N_DN=nrow(subset(sig,log2FoldChange<0))
DET=nrow(de)
HEADER=paste(N_SIG,"@5%FDR,", N_UP, "up", N_DN, "dn", DET, "detected")
plot(de$log2FoldChange,-log10(de$pval),cex=0.5,pch=19,col="darkgray",
main=name, xlab="log2 FC", ylab="-log10 pval")
mtext(HEADER)
grid()
points(sig$log2FoldChange,-log10(sig$pval),cex=0.5,pch=19,col="red")
}
make_heatmap <- function(de,name,myss,mx,n=30){
colfunc <- colorRampPalette(c("blue", "white", "red"))
mxn <- mx/rowSums(mx)*1000000
x <- mxn[which(rownames(mxn) %in% rownames(head(de,n))),]
heatmap.2(as.matrix(x),trace="none",col=colfunc(25),scale="row", margins = c(10,15), cexRow=0.7,
main=paste("Top ranked",n,"genes in",name) , ColSideColors = myss$cols )
}
maplot(dge1,"Cont1: Effect of KO in untr cells")
make_volcano(dge1,"Cont1: Effect of KO in untr cells")
make_heatmap(dge1,"Cont1: Effect of KO in untr cells",ss1,xx1,n=50)
maplot(dge2,"Cont2: Effect of KO in trt cells")
make_volcano(dge2,"Cont2: Effect of KO in trt cells")
make_heatmap(dge2,"Cont2: Effect of KO in trt cells",ss2,xx2,n=50)
maplot(dge3,"Cont3: Effect of drug in WT cells")
make_volcano(dge3,"Cont3: Effect of drug in WT cells")
make_heatmap(dge3,"Cont3: Effect of drug in WT cells",ss3,xx3,n=50)
maplot(dge4,"Cont4: Effect of drug in KO cells")
make_volcano(dge4,"Cont4: Effect of drug in KO cells")
make_heatmap(dge4,"Cont4: Effect of drug in KO cells",ss4,xx4,n=50)
Euler diagrams of DEGs
All DEGs
par(cex.main=0.5)
par(mar=c(2,2,2,2))
v0 <- list("C1 up"=d1up,"C1 dn"=d1dn, "C2 up" = d2up, "C2 dn"=d2dn, "C3 up"=d3up, "C3 dn"=d3dn, "C4 up"=d4up, "C4 dn"=d4dn)
plot(euler(v0),quantities = TRUE, edges = "gray", main="DEGs")
Now with just the DEGs from contrasts 1 and 4.
par(cex.main=0.5)
par(mar=c(2,2,2,2))
v0 <- list("KO up"=d1up,"KO dn"=d1dn, "drug up"=d4up, "drug dn"=d4dn)
plot(euler(v0),quantities = TRUE, edges = "gray", main="DEGs")
Single contrast pathway analysis with mitch
Firstly need to conduct mitch enrichment analysis for each contrast separately.
#download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip", destfile="ReactomePathways.gmt.zip")
#unzip("ReactomePathways.gmt.zip")
genesets <- gmt_import("ReactomePathways.gmt")
# gene table
gt <- as.data.frame(rownames(xx))
gt$gene <- sapply(strsplit(gt[,1]," "),"[[",2)
m1 <- mitch_import(dge1, 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 = 25258## Note: no. genes in output = 25203## Note: estimated proportion of input genes in output = 0.998mres1 <- mitch_calc(m1, genesets, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(mres1$enrichment_result,20) %>% kbl(caption = "Top gene pathway differences in contrast 1") %>% kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 1008 | Receptor Mediated Mitophagy | 11 | 0.0000080 | -0.7772236 | 0.0000639 |
| 720 | Mucopolysaccharidoses | 11 | 0.0000117 | 0.7629839 | 0.0000908 |
| 670 | Membrane binding and targetting of GAG proteins | 14 | 0.0000402 | -0.6337233 | 0.0002665 |
| 1264 | Synthesis And Processing Of GAG, GAGPOL Polyproteins | 14 | 0.0000402 | -0.6337233 | 0.0002665 |
| 1411 | Unwinding of DNA | 12 | 0.0002829 | 0.6051698 | 0.0014768 |
| 93 | Assembly Of The HIV Virion | 16 | 0.0000324 | -0.6000020 | 0.0002240 |
| 121 | Biotin transport and metabolism | 11 | 0.0011127 | 0.5676477 | 0.0047067 |
| 1276 | Synthesis of active ubiquitin: roles of E1 and E2 enzymes | 30 | 0.0000001 | -0.5630345 | 0.0000016 |
| 1115 | SLBP Dependent Processing of Replication-Dependent Histone Pre-mRNAs | 11 | 0.0018488 | -0.5420841 | 0.0073752 |
| 1166 | Signaling by Activin | 15 | 0.0003295 | -0.5354190 | 0.0016612 |
| 710 | Mitophagy | 28 | 0.0000011 | -0.5309604 | 0.0000116 |
| 1086 | Response of EIF2AK1 (HRI) to heme deficiency | 15 | 0.0005010 | -0.5189191 | 0.0023789 |
| 55 | Activation of the mRNA upon binding of the cap-binding complex and eIFs, and subsequent binding to 43S | 58 | 0.0000000 | -0.5187746 | 0.0000000 |
| 538 | Hyaluronan uptake and degradation | 10 | 0.0047720 | 0.5153177 | 0.0167247 |
| 572 | Insertion of tail-anchored proteins into the endoplasmic reticulum membrane | 22 | 0.0000312 | -0.5127387 | 0.0002190 |
| 1376 | Translation initiation complex formation | 57 | 0.0000000 | -0.5104645 | 0.0000000 |
| 644 | Lysosphingolipid and LPA receptors | 12 | 0.0024056 | -0.5059280 | 0.0092214 |
| 1462 | rRNA modification in the nucleus and cytosol | 60 | 0.0000000 | -0.5021477 | 0.0000000 |
| 517 | HIV elongation arrest and recovery | 32 | 0.0000009 | -0.5012440 | 0.0000096 |
| 848 | Pausing and recovery of HIV elongation | 32 | 0.0000009 | -0.5012440 | 0.0000096 |
m1top <- subset(mres1$enrichment_result,p.adjustANOVA<0.05)
m1up <- subset(m1top,s.dist>0)$set
m1dn <- subset(m1top,s.dist<0)$set
mitch_report(mres1,outfile="mitch1.html",overwrite=TRUE)## Note: overwriting existing report## Dataset saved as " /tmp/RtmpMZfPOC/mitch1.rds ".##
##
## processing file: mitch.Rmd##
|
| | 0%
|
|.. | 3%
## inline R code fragments
##
##
|
|.... | 6%
## label: checklibraries (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|...... | 9%
## ordinary text without R code
##
##
|
|........ | 12%
## label: peek (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.......... | 15%
## ordinary text without R code
##
##
|
|............ | 18%
## label: metrics (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.............. | 21%
## ordinary text without R code
##
##
|
|................ | 24%
## label: scatterplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ message : logi FALSE
## $ warning : logi FALSE##
|
|................... | 26%
## ordinary text without R code
##
##
|
|..................... | 29%
## label: contourplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ warning : logi FALSE
## $ message : logi FALSE## Contour plot does not apply to unidimensional analysis.##
|
|....................... | 32%
## ordinary text without R code
##
##
|
|......................... | 35%
## label: input_geneset_metrics1 (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|........................... | 38%
## ordinary text without R code
##
##
|
|............................. | 41%
## label: input_geneset_metrics2 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"##
|
|............................... | 44%
## ordinary text without R code
##
##
|
|................................. | 47%
## label: input_geneset_metrics3 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ message : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7##
|
|................................... | 50%
## ordinary text without R code
##
##
|
|..................................... | 53%
## label: echart1d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|....................................... | 56%
## label: echart2d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|......................................... | 59%
## ordinary text without R code
##
##
|
|........................................... | 62%
## label: heatmap (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 10
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|............................................. | 65%
## ordinary text without R code
##
##
|
|............................................... | 68%
## label: effectsize (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|................................................. | 71%
## ordinary text without R code
##
##
|
|................................................... | 74%
## label: results_table (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|...................................................... | 76%
## ordinary text without R code
##
##
|
|........................................................ | 79%
## label: results_table_complete (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|.......................................................... | 82%
## ordinary text without R code
##
##
|
|............................................................ | 85%
## label: detailed_geneset_reports1d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE##
|
|.............................................................. | 88%
## ordinary text without R code
##
##
|
|................................................................ | 91%
## label: detailed_geneset_reports2d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 5
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE
##
##
|
|.................................................................. | 94%
## ordinary text without R code
##
##
|
|.................................................................... | 97%
## label: session_info (with options)
## List of 3
## $ include: logi TRUE
## $ echo : logi TRUE
## $ results: chr "markup"
##
##
|
|......................................................................| 100%
## ordinary text without R code## output file: /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md## /home/mdz/anaconda3/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpMZfPOC/mitch_report.html --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --standalone --section-divs --template /usr/lib/R/library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --include-in-header /tmp/RtmpMZfPOC/rmarkdown-str179303fa111de.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'##
## Output created: /tmp/RtmpMZfPOC/mitch_report.html## [1] TRUEwrite.table(mres1$enrichment_result,file="mitch1.tsv",quote=FALSE,sep="\t",row.names=FALSE)
m2 <- mitch_import(dge2, 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 = 27340## Note: no. genes in output = 27278## Note: estimated proportion of input genes in output = 0.998mres2 <- mitch_calc(m2, genesets, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(mres2$enrichment_result,20) %>% kbl(caption = "Top gene pathway differences in contrast 2") %>% kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 1427 | Unwinding of DNA | 12 | 0.0000101 | 0.7359532 | 0.0000679 |
| 640 | Leading Strand Synthesis | 14 | 0.0000029 | 0.7219567 | 0.0000217 |
| 893 | Polymerase switching | 14 | 0.0000029 | 0.7219567 | 0.0000217 |
| 1102 | Response to metal ions | 10 | 0.0000778 | 0.7214390 | 0.0004521 |
| 458 | G1/S-Specific Transcription | 29 | 0.0000000 | 0.6885884 | 0.0000000 |
| 216 | Condensation of Prometaphase Chromosomes | 11 | 0.0000792 | 0.6871343 | 0.0004567 |
| 1345 | Tetrahydrobiopterin (BH4) synthesis, recycling, salvage and regulation | 10 | 0.0002313 | 0.6723265 | 0.0011212 |
| 423 | Formation of ATP by chemiosmotic coupling | 18 | 0.0000009 | 0.6681544 | 0.0000075 |
| 1130 | SLBP independent Processing of Histone Pre-mRNAs | 10 | 0.0003221 | 0.6567625 | 0.0015025 |
| 1129 | SLBP Dependent Processing of Replication-Dependent Histone Pre-mRNAs | 11 | 0.0001693 | 0.6547542 | 0.0008750 |
| 260 | DNA strand elongation | 32 | 0.0000000 | 0.6539698 | 0.0000000 |
| 1397 | Translesion synthesis by REV1 | 16 | 0.0000094 | 0.6395441 | 0.0000640 |
| 1396 | Translesion synthesis by POLK | 17 | 0.0000057 | 0.6353722 | 0.0000407 |
| 422 | Folding of actin by CCT/TriC | 10 | 0.0005533 | 0.6306000 | 0.0024504 |
| 1118 | SCF(Skp2)-mediated degradation of p27/p21 | 60 | 0.0000000 | 0.6273814 | 0.0000000 |
| 1351 | The role of GTSE1 in G2/M progression after G2 checkpoint | 58 | 0.0000000 | 0.6165776 | 0.0000000 |
| 686 | Metabolism of cofactors | 19 | 0.0000036 | 0.6139295 | 0.0000265 |
| 1050 | Regulation of PTEN mRNA translation | 12 | 0.0002439 | -0.6115125 | 0.0011746 |
| 1424 | Ubiquitin-dependent degradation of Cyclin D | 50 | 0.0000000 | 0.6083723 | 0.0000000 |
| 636 | Lagging Strand Synthesis | 20 | 0.0000029 | 0.6041896 | 0.0000217 |
m2top <- subset(mres2$enrichment_result,p.adjustANOVA<0.05)
m2up <- subset(m2top,s.dist>0)$set
m2dn <- subset(m2top,s.dist<0)$set
mitch_report(mres2,outfile="mitch2.html",overwrite=TRUE)## Note: overwriting existing report## Dataset saved as " /tmp/RtmpMZfPOC/mitch2.rds ".##
##
## processing file: mitch.Rmd##
|
| | 0%
|
|.. | 3%
## inline R code fragments
##
##
|
|.... | 6%
## label: checklibraries (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|...... | 9%
## ordinary text without R code
##
##
|
|........ | 12%
## label: peek (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.......... | 15%
## ordinary text without R code
##
##
|
|............ | 18%
## label: metrics (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.............. | 21%
## ordinary text without R code
##
##
|
|................ | 24%
## label: scatterplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ message : logi FALSE
## $ warning : logi FALSE##
|
|................... | 26%
## ordinary text without R code
##
##
|
|..................... | 29%
## label: contourplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ warning : logi FALSE
## $ message : logi FALSE## Contour plot does not apply to unidimensional analysis.##
|
|....................... | 32%
## ordinary text without R code
##
##
|
|......................... | 35%
## label: input_geneset_metrics1 (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|........................... | 38%
## ordinary text without R code
##
##
|
|............................. | 41%
## label: input_geneset_metrics2 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"##
|
|............................... | 44%
## ordinary text without R code
##
##
|
|................................. | 47%
## label: input_geneset_metrics3 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ message : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7##
|
|................................... | 50%
## ordinary text without R code
##
##
|
|..................................... | 53%
## label: echart1d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|....................................... | 56%
## label: echart2d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|......................................... | 59%
## ordinary text without R code
##
##
|
|........................................... | 62%
## label: heatmap (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 10
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|............................................. | 65%
## ordinary text without R code
##
##
|
|............................................... | 68%
## label: effectsize (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|................................................. | 71%
## ordinary text without R code
##
##
|
|................................................... | 74%
## label: results_table (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|...................................................... | 76%
## ordinary text without R code
##
##
|
|........................................................ | 79%
## label: results_table_complete (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|.......................................................... | 82%
## ordinary text without R code
##
##
|
|............................................................ | 85%
## label: detailed_geneset_reports1d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE##
|
|.............................................................. | 88%
## ordinary text without R code
##
##
|
|................................................................ | 91%
## label: detailed_geneset_reports2d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 5
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE
##
##
|
|.................................................................. | 94%
## ordinary text without R code
##
##
|
|.................................................................... | 97%
## label: session_info (with options)
## List of 3
## $ include: logi TRUE
## $ echo : logi TRUE
## $ results: chr "markup"
##
##
|
|......................................................................| 100%
## ordinary text without R code## output file: /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md## /home/mdz/anaconda3/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpMZfPOC/mitch_report.html --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --standalone --section-divs --template /usr/lib/R/library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --include-in-header /tmp/RtmpMZfPOC/rmarkdown-str1793022e565d.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'##
## Output created: /tmp/RtmpMZfPOC/mitch_report.html## [1] TRUEwrite.table(mres2$enrichment_result,file="mitch2.tsv",quote=FALSE,sep="\t",row.names=FALSE)
m3 <- mitch_import(dge3, 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 = 26801## Note: no. genes in output = 26741## Note: estimated proportion of input genes in output = 0.998mres3 <- mitch_calc(m3, genesets, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(mres3$enrichment_result,20) %>% kbl(caption = "Top gene pathway differences in contrast 3") %>% kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 137 | CLEC7A (Dectin-1) induces NFAT activation | 11 | 1.70e-06 | -0.8325273 | 0.0000081 |
| 375 | Establishment of Sister Chromatid Cohesion | 11 | 2.30e-06 | -0.8231201 | 0.0000102 |
| 207 | Cohesin Loading onto Chromatin | 10 | 7.30e-06 | -0.8190042 | 0.0000294 |
| 360 | ERKs are inactivated | 13 | 4.00e-07 | -0.8096839 | 0.0000023 |
| 1108 | S33 mutants of beta-catenin aren’t phosphorylated | 15 | 1.00e-07 | -0.8083414 | 0.0000004 |
| 1109 | S37 mutants of beta-catenin aren’t phosphorylated | 15 | 1.00e-07 | -0.8083414 | 0.0000004 |
| 1110 | S45 mutants of beta-catenin aren’t phosphorylated | 15 | 1.00e-07 | -0.8083414 | 0.0000004 |
| 1180 | Signaling by CTNNB1 phospho-site mutants | 15 | 1.00e-07 | -0.8083414 | 0.0000004 |
| 1207 | Signaling by GSK3beta mutants | 15 | 1.00e-07 | -0.8083414 | 0.0000004 |
| 1298 | T41 mutants of beta-catenin aren’t phosphorylated | 15 | 1.00e-07 | -0.8083414 | 0.0000004 |
| 506 | Golgi Cisternae Pericentriolar Stack Reorganization | 14 | 2.00e-07 | -0.8059533 | 0.0000011 |
| 11 | APC truncation mutants have impaired AXIN binding | 14 | 2.00e-07 | -0.7983954 | 0.0000013 |
| 24 | AXIN missense mutants destabilize the destruction complex | 14 | 2.00e-07 | -0.7983954 | 0.0000013 |
| 1173 | Signaling by AMER1 mutants | 14 | 2.00e-07 | -0.7983954 | 0.0000013 |
| 1174 | Signaling by APC mutants | 14 | 2.00e-07 | -0.7983954 | 0.0000013 |
| 1175 | Signaling by AXIN mutants | 14 | 2.00e-07 | -0.7983954 | 0.0000013 |
| 1414 | Truncations of AMER1 destabilize the destruction complex | 14 | 2.00e-07 | -0.7983954 | 0.0000013 |
| 1209 | Signaling by Hippo | 20 | 0.00e+00 | -0.7935032 | 0.0000000 |
| 664 | MET activates RAP1 and RAC1 | 11 | 5.40e-06 | -0.7917695 | 0.0000227 |
| 659 | MAPK3 (ERK1) activation | 10 | 2.78e-05 | -0.7651341 | 0.0001004 |
m3top <- subset(mres3$enrichment_result,p.adjustANOVA<0.05)
m3up <- subset(m3top,s.dist>0)$set
m3dn <- subset(m3top,s.dist<0)$set
mitch_report(mres3,outfile="mitch3.html",overwrite=TRUE)## Note: overwriting existing report## Dataset saved as " /tmp/RtmpMZfPOC/mitch3.rds ".##
##
## processing file: mitch.Rmd##
|
| | 0%
|
|.. | 3%
## inline R code fragments
##
##
|
|.... | 6%
## label: checklibraries (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|...... | 9%
## ordinary text without R code
##
##
|
|........ | 12%
## label: peek (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.......... | 15%
## ordinary text without R code
##
##
|
|............ | 18%
## label: metrics (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.............. | 21%
## ordinary text without R code
##
##
|
|................ | 24%
## label: scatterplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ message : logi FALSE
## $ warning : logi FALSE##
|
|................... | 26%
## ordinary text without R code
##
##
|
|..................... | 29%
## label: contourplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ warning : logi FALSE
## $ message : logi FALSE## Contour plot does not apply to unidimensional analysis.##
|
|....................... | 32%
## ordinary text without R code
##
##
|
|......................... | 35%
## label: input_geneset_metrics1 (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|........................... | 38%
## ordinary text without R code
##
##
|
|............................. | 41%
## label: input_geneset_metrics2 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"##
|
|............................... | 44%
## ordinary text without R code
##
##
|
|................................. | 47%
## label: input_geneset_metrics3 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ message : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7##
|
|................................... | 50%
## ordinary text without R code
##
##
|
|..................................... | 53%
## label: echart1d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|....................................... | 56%
## label: echart2d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|......................................... | 59%
## ordinary text without R code
##
##
|
|........................................... | 62%
## label: heatmap (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 10
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|............................................. | 65%
## ordinary text without R code
##
##
|
|............................................... | 68%
## label: effectsize (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|................................................. | 71%
## ordinary text without R code
##
##
|
|................................................... | 74%
## label: results_table (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|...................................................... | 76%
## ordinary text without R code
##
##
|
|........................................................ | 79%
## label: results_table_complete (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|.......................................................... | 82%
## ordinary text without R code
##
##
|
|............................................................ | 85%
## label: detailed_geneset_reports1d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE##
|
|.............................................................. | 88%
## ordinary text without R code
##
##
|
|................................................................ | 91%
## label: detailed_geneset_reports2d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 5
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE
##
##
|
|.................................................................. | 94%
## ordinary text without R code
##
##
|
|.................................................................... | 97%
## label: session_info (with options)
## List of 3
## $ include: logi TRUE
## $ echo : logi TRUE
## $ results: chr "markup"
##
##
|
|......................................................................| 100%
## ordinary text without R code## output file: /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md## /home/mdz/anaconda3/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpMZfPOC/mitch_report.html --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --standalone --section-divs --template /usr/lib/R/library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --include-in-header /tmp/RtmpMZfPOC/rmarkdown-str17930cf8a15e.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'##
## Output created: /tmp/RtmpMZfPOC/mitch_report.html## [1] TRUEwrite.table(mres3$enrichment_result,file="mitch3.tsv",quote=FALSE,sep="\t",row.names=FALSE)
m4 <- mitch_import(dge4, 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 = 26222## Note: no. genes in output = 26163## Note: estimated proportion of input genes in output = 0.998mres4 <- mitch_calc(m4, genesets, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(mres4$enrichment_result,20) %>% kbl(caption = "Top gene pathway differences in contrast 4") %>% kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 134 | CLEC7A (Dectin-1) induces NFAT activation | 11 | 1.0e-06 | -0.8526586 | 6.10e-06 |
| 203 | Cohesin Loading onto Chromatin | 10 | 3.0e-06 | -0.8523305 | 1.63e-05 |
| 630 | Laminin interactions | 23 | 0.0e+00 | -0.7933735 | 0.00e+00 |
| 1103 | S33 mutants of beta-catenin aren’t phosphorylated | 15 | 2.0e-07 | -0.7825965 | 1.20e-06 |
| 1104 | S37 mutants of beta-catenin aren’t phosphorylated | 15 | 2.0e-07 | -0.7825965 | 1.20e-06 |
| 1105 | S45 mutants of beta-catenin aren’t phosphorylated | 15 | 2.0e-07 | -0.7825965 | 1.20e-06 |
| 1175 | Signaling by CTNNB1 phospho-site mutants | 15 | 2.0e-07 | -0.7825965 | 1.20e-06 |
| 1202 | Signaling by GSK3beta mutants | 15 | 2.0e-07 | -0.7825965 | 1.20e-06 |
| 1293 | T41 mutants of beta-catenin aren’t phosphorylated | 15 | 2.0e-07 | -0.7825965 | 1.20e-06 |
| 1204 | Signaling by Hippo | 20 | 0.0e+00 | -0.7785564 | 0.00e+00 |
| 1268 | Syndecan interactions | 20 | 0.0e+00 | -0.7770990 | 0.00e+00 |
| 658 | MET activates RAP1 and RAC1 | 11 | 8.9e-06 | -0.7734156 | 4.31e-05 |
| 11 | APC truncation mutants have impaired AXIN binding | 14 | 6.0e-07 | -0.7679398 | 4.20e-06 |
| 24 | AXIN missense mutants destabilize the destruction complex | 14 | 6.0e-07 | -0.7679398 | 4.20e-06 |
| 1168 | Signaling by AMER1 mutants | 14 | 6.0e-07 | -0.7679398 | 4.20e-06 |
| 1169 | Signaling by APC mutants | 14 | 6.0e-07 | -0.7679398 | 4.20e-06 |
| 1170 | Signaling by AXIN mutants | 14 | 6.0e-07 | -0.7679398 | 4.20e-06 |
| 1410 | Truncations of AMER1 destabilize the destruction complex | 14 | 6.0e-07 | -0.7679398 | 4.20e-06 |
| 357 | ERKs are inactivated | 13 | 1.7e-06 | -0.7668245 | 9.70e-06 |
| 657 | MET activates PTK2 signaling | 18 | 0.0e+00 | -0.7616795 | 2.00e-07 |
m4top <- subset(mres4$enrichment_result,p.adjustANOVA<0.05)
m4up <- subset(m4top,s.dist>0)$set
m4dn <- subset(m4top,s.dist<0)$set
mitch_report(mres4,outfile="mitch4.html",overwrite=TRUE)## Note: overwriting existing report## Dataset saved as " /tmp/RtmpMZfPOC/mitch4.rds ".##
##
## processing file: mitch.Rmd##
|
| | 0%
|
|.. | 3%
## inline R code fragments
##
##
|
|.... | 6%
## label: checklibraries (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|...... | 9%
## ordinary text without R code
##
##
|
|........ | 12%
## label: peek (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.......... | 15%
## ordinary text without R code
##
##
|
|............ | 18%
## label: metrics (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.............. | 21%
## ordinary text without R code
##
##
|
|................ | 24%
## label: scatterplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ message : logi FALSE
## $ warning : logi FALSE##
|
|................... | 26%
## ordinary text without R code
##
##
|
|..................... | 29%
## label: contourplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ warning : logi FALSE
## $ message : logi FALSE## Contour plot does not apply to unidimensional analysis.##
|
|....................... | 32%
## ordinary text without R code
##
##
|
|......................... | 35%
## label: input_geneset_metrics1 (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|........................... | 38%
## ordinary text without R code
##
##
|
|............................. | 41%
## label: input_geneset_metrics2 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"##
|
|............................... | 44%
## ordinary text without R code
##
##
|
|................................. | 47%
## label: input_geneset_metrics3 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ message : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7##
|
|................................... | 50%
## ordinary text without R code
##
##
|
|..................................... | 53%
## label: echart1d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|....................................... | 56%
## label: echart2d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|......................................... | 59%
## ordinary text without R code
##
##
|
|........................................... | 62%
## label: heatmap (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 10
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|............................................. | 65%
## ordinary text without R code
##
##
|
|............................................... | 68%
## label: effectsize (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|................................................. | 71%
## ordinary text without R code
##
##
|
|................................................... | 74%
## label: results_table (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|...................................................... | 76%
## ordinary text without R code
##
##
|
|........................................................ | 79%
## label: results_table_complete (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|.......................................................... | 82%
## ordinary text without R code
##
##
|
|............................................................ | 85%
## label: detailed_geneset_reports1d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE##
|
|.............................................................. | 88%
## ordinary text without R code
##
##
|
|................................................................ | 91%
## label: detailed_geneset_reports2d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 5
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE
##
##
|
|.................................................................. | 94%
## ordinary text without R code
##
##
|
|.................................................................... | 97%
## label: session_info (with options)
## List of 3
## $ include: logi TRUE
## $ echo : logi TRUE
## $ results: chr "markup"
##
##
|
|......................................................................| 100%
## ordinary text without R code## output file: /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md## /home/mdz/anaconda3/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpMZfPOC/mitch_report.html --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --standalone --section-divs --template /usr/lib/R/library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --include-in-header /tmp/RtmpMZfPOC/rmarkdown-str179307b80719f.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'##
## Output created: /tmp/RtmpMZfPOC/mitch_report.html## [1] TRUEwrite.table(mres4$enrichment_result,file="mitch4.tsv",quote=FALSE,sep="\t",row.names=FALSE)Now an Euler diagram of differentially regulated gene sets from contrasts 1 and 4.
par(cex.main=0.5)
par(mar=c(2,2,2,2))
v0 <- list("KO up"=m1up,"KO dn"=m1dn, "drug up"=m4up, "drug dn"=m4dn)
plot(euler(v0),quantities = TRUE, edges = "gray", main="Reactome pathways")
Multi contrast pathway analysis with mitch
I’m going to focus on the two contrasts which I think are the most important WRT mito disease; #1 and #4. The effect of the KO, and the effect of the drug when the gene is knocked out.
mm <- list("KO"=dge1,"drug"=dge4)
m <- mitch_import(mm, DEtype="deseq2",geneTable=gt)## Note: Mean no. genes in input = 25740## Note: no. genes in output = 24473## Note: estimated proportion of input genes in output = 0.951mres <- mitch_calc(m, genesets, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(mres$enrichment_result,20) %>%
kbl(caption = "Top gene pathway differences in contrast 1 and 4") %>% kable_paper("hover", full_width = F)| set | setSize | pMANOVA | s.KO | s.drug | p.KO | p.drug | s.dist | SD | p.adjustMANOVA | |
|---|---|---|---|---|---|---|---|---|---|---|
| 133 | CLEC7A (Dectin-1) induces NFAT activation | 11 | 3.40e-06 | -0.3079507 | -0.8424792 | 0.0769588 | 0.0000013 | 0.8969977 | 0.3779687 | 1.31e-05 |
| 202 | Cohesin Loading onto Chromatin | 10 | 1.94e-05 | -0.1962392 | -0.8421289 | 0.2825607 | 0.0000040 | 0.8646913 | 0.4567130 | 6.18e-05 |
| 650 | MET activates RAP1 and RAC1 | 11 | 1.26e-05 | -0.3986517 | -0.7577616 | 0.0220462 | 0.0000135 | 0.8562277 | 0.2539290 | 4.27e-05 |
| 1256 | Syndecan interactions | 19 | 0.00e+00 | 0.2156918 | -0.8161446 | 0.1035746 | 0.0000000 | 0.8441652 | 0.7296185 | 0.00e+00 |
| 844 | Peptide chain elongation | 87 | 0.00e+00 | -0.4224866 | 0.7146629 | 0.0000000 | 0.0000000 | 0.8302036 | 0.8040862 | 0.00e+00 |
| 1415 | Viral mRNA Translation | 87 | 0.00e+00 | -0.4204815 | 0.7154708 | 0.0000000 | 0.0000000 | 0.8298814 | 0.8032396 | 0.00e+00 |
| 11 | APC truncation mutants have impaired AXIN binding | 14 | 1.30e-06 | -0.3476780 | -0.7519055 | 0.0242864 | 0.0000011 | 0.8283972 | 0.2858320 | 5.60e-06 |
| 24 | AXIN missense mutants destabilize the destruction complex | 14 | 1.30e-06 | -0.3476780 | -0.7519055 | 0.0242864 | 0.0000011 | 0.8283972 | 0.2858320 | 5.60e-06 |
| 1156 | Signaling by AMER1 mutants | 14 | 1.30e-06 | -0.3476780 | -0.7519055 | 0.0242864 | 0.0000011 | 0.8283972 | 0.2858320 | 5.60e-06 |
| 1157 | Signaling by APC mutants | 14 | 1.30e-06 | -0.3476780 | -0.7519055 | 0.0242864 | 0.0000011 | 0.8283972 | 0.2858320 | 5.60e-06 |
| 1158 | Signaling by AXIN mutants | 14 | 1.30e-06 | -0.3476780 | -0.7519055 | 0.0242864 | 0.0000011 | 0.8283972 | 0.2858320 | 5.60e-06 |
| 1396 | Truncations of AMER1 destabilize the destruction complex | 14 | 1.30e-06 | -0.3476780 | -0.7519055 | 0.0242864 | 0.0000011 | 0.8283972 | 0.2858320 | 5.60e-06 |
| 622 | Laminin interactions | 23 | 0.00e+00 | 0.2685267 | -0.7821464 | 0.0257806 | 0.0000000 | 0.8269580 | 0.7429381 | 0.00e+00 |
| 351 | ERKs are inactivated | 13 | 4.30e-06 | -0.3333229 | -0.7507139 | 0.0374309 | 0.0000028 | 0.8213863 | 0.2951400 | 1.60e-05 |
| 1001 | Receptor Mediated Mitophagy | 11 | 2.63e-05 | -0.7759568 | -0.2626486 | 0.0000083 | 0.1314471 | 0.8192028 | 0.3629637 | 8.10e-05 |
| 1091 | S33 mutants of beta-catenin aren’t phosphorylated | 15 | 8.00e-07 | -0.2603974 | -0.7675743 | 0.0807712 | 0.0000003 | 0.8105413 | 0.3586282 | 3.50e-06 |
| 1092 | S37 mutants of beta-catenin aren’t phosphorylated | 15 | 8.00e-07 | -0.2603974 | -0.7675743 | 0.0807712 | 0.0000003 | 0.8105413 | 0.3586282 | 3.50e-06 |
| 1093 | S45 mutants of beta-catenin aren’t phosphorylated | 15 | 8.00e-07 | -0.2603974 | -0.7675743 | 0.0807712 | 0.0000003 | 0.8105413 | 0.3586282 | 3.50e-06 |
| 1163 | Signaling by CTNNB1 phospho-site mutants | 15 | 8.00e-07 | -0.2603974 | -0.7675743 | 0.0807712 | 0.0000003 | 0.8105413 | 0.3586282 | 3.50e-06 |
| 1190 | Signaling by GSK3beta mutants | 15 | 8.00e-07 | -0.2603974 | -0.7675743 | 0.0807712 | 0.0000003 | 0.8105413 | 0.3586282 | 3.50e-06 |
mitch_report(mres,outfile="mitchmulti_effect.html",overwrite=TRUE)## Note: overwriting existing report## Dataset saved as " /tmp/RtmpMZfPOC/mitchmulti_effect.rds ".##
##
## processing file: mitch.Rmd##
|
| | 0%
|
|.. | 3%
## inline R code fragments
##
##
|
|.... | 6%
## label: checklibraries (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|...... | 9%
## ordinary text without R code
##
##
|
|........ | 12%
## label: peek (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.......... | 15%
## ordinary text without R code
##
##
|
|............ | 18%
## label: metrics (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.............. | 21%
## ordinary text without R code
##
##
|
|................ | 24%
## label: scatterplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ message : logi FALSE
## $ warning : logi FALSE##
|
|................... | 26%
## ordinary text without R code
##
##
|
|..................... | 29%
## label: contourplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ warning : logi FALSE
## $ message : logi FALSE##
|
|....................... | 32%
## ordinary text without R code
##
##
|
|......................... | 35%
## label: input_geneset_metrics1 (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|........................... | 38%
## ordinary text without R code
##
##
|
|............................. | 41%
## label: input_geneset_metrics2 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"##
|
|............................... | 44%
## ordinary text without R code
##
##
|
|................................. | 47%
## label: input_geneset_metrics3 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ message : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7##
|
|................................... | 50%
## ordinary text without R code
##
##
|
|..................................... | 53%
## label: echart1d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|....................................... | 56%
## label: echart2d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|......................................... | 59%
## ordinary text without R code
##
##
|
|........................................... | 62%
## label: heatmap (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 10
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE##
|
|............................................. | 65%
## ordinary text without R code
##
##
|
|............................................... | 68%
## label: effectsize (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|................................................. | 71%
## ordinary text without R code
##
##
|
|................................................... | 74%
## label: results_table (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|...................................................... | 76%
## ordinary text without R code
##
##
|
|........................................................ | 79%
## label: results_table_complete (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|.......................................................... | 82%
## ordinary text without R code
##
##
|
|............................................................ | 85%
## label: detailed_geneset_reports1d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE
##
##
|
|.............................................................. | 88%
## ordinary text without R code
##
##
|
|................................................................ | 91%
## label: detailed_geneset_reports2d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 5
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE##
|
|.................................................................. | 94%
## ordinary text without R code
##
##
|
|.................................................................... | 97%
## label: session_info (with options)
## List of 3
## $ include: logi TRUE
## $ echo : logi TRUE
## $ results: chr "markup"
##
##
|
|......................................................................| 100%
## ordinary text without R code## output file: /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md## /home/mdz/anaconda3/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpMZfPOC/mitch_report.html --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --standalone --section-divs --template /usr/lib/R/library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --include-in-header /tmp/RtmpMZfPOC/rmarkdown-str1793076bd7e3.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'##
## Output created: /tmp/RtmpMZfPOC/mitch_report.html## [1] TRUEwrite.table(mres$enrichment_result,file="mitch_multi_effect.tsv",quote=FALSE,sep="\t",row.names=FALSE)
n=40
mx <- as.matrix(mres$enrichment_result[1:n,4:5])
rownames(mx) <- mres$set
rownames(mx) <- mres$enrichment_result$set[1:n]
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(mx,trace="none",col=colfunc(25),scale="none", margins = c(5,15), cexRow=0.7,
cexCol=0.9, main=paste("Top ranked Reactomes") )
Mitch can also be used to search for discordant pathways. Gene sets are prioritised by standard deviation.
mm <- list("KO"=dge1,"drug"=dge4)
m <- mitch_import(mm, DEtype="deseq2",geneTable=gt)## Note: Mean no. genes in input = 25740## Note: no. genes in output = 24473## Note: estimated proportion of input genes in output = 0.951mres <- mitch_calc(m, genesets, priority="SD")## Note: Prioritisation by SD after selecting sets with
## p.adjustMANOVA<=0.05.head(mres$enrichment_result,20) %>% kbl(caption = "Top gene pathway differences in contrast 1 and 4") %>% kable_paper("hover", full_width = F)| set | setSize | pMANOVA | s.KO | s.drug | p.KO | p.drug | s.dist | SD | p.adjustMANOVA | |
|---|---|---|---|---|---|---|---|---|---|---|
| 844 | Peptide chain elongation | 87 | 0.00e+00 | -0.4224866 | 0.7146629 | 0.0000000 | 0.0000000 | 0.8302036 | 0.8040862 | 0.00e+00 |
| 1415 | Viral mRNA Translation | 87 | 0.00e+00 | -0.4204815 | 0.7154708 | 0.0000000 | 0.0000000 | 0.8298814 | 0.8032396 | 0.00e+00 |
| 419 | Formation of a pool of free 40S subunits | 99 | 0.00e+00 | -0.4444643 | 0.6516333 | 0.0000000 | 0.0000000 | 0.7887804 | 0.7750580 | 0.00e+00 |
| 371 | Eukaryotic Translation Termination | 91 | 0.00e+00 | -0.4020927 | 0.6842996 | 0.0000000 | 0.0000000 | 0.7936904 | 0.7681954 | 0.00e+00 |
| 369 | Eukaryotic Translation Elongation | 92 | 0.00e+00 | -0.4120627 | 0.6735448 | 0.0000000 | 0.0000000 | 0.7895937 | 0.7676404 | 0.00e+00 |
| 268 | Defective EXT1 causes exostoses 1, TRPS2 and CHDS | 13 | 1.40e-06 | 0.3576640 | -0.7141141 | 0.0255478 | 0.0000082 | 0.7986755 | 0.7578616 | 5.90e-06 |
| 269 | Defective EXT2 causes exostoses 2 | 13 | 1.40e-06 | 0.3576640 | -0.7141141 | 0.0255478 | 0.0000082 | 0.7986755 | 0.7578616 | 5.90e-06 |
| 622 | Laminin interactions | 23 | 0.00e+00 | 0.2685267 | -0.7821464 | 0.0257806 | 0.0000000 | 0.8269580 | 0.7429381 | 0.00e+00 |
| 1080 | Response of EIF2AK4 (GCN2) to amino acid deficiency | 99 | 0.00e+00 | -0.4230307 | 0.6251677 | 0.0000000 | 0.0000000 | 0.7548441 | 0.7411882 | 0.00e+00 |
| 777 | Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC) | 93 | 0.00e+00 | -0.4238720 | 0.6173631 | 0.0000000 | 0.0000000 | 0.7488689 | 0.7362644 | 0.00e+00 |
| 1134 | Selenocysteine synthesis | 91 | 0.00e+00 | -0.3552729 | 0.6805615 | 0.0000000 | 0.0000000 | 0.7677127 | 0.7324455 | 0.00e+00 |
| 425 | Formation of the ternary complex, and subsequently, the 43S complex | 50 | 0.00e+00 | -0.5007952 | 0.5323556 | 0.0000000 | 0.0000000 | 0.7308887 | 0.7305479 | 0.00e+00 |
| 1256 | Syndecan interactions | 19 | 0.00e+00 | 0.2156918 | -0.8161446 | 0.1035746 | 0.0000000 | 0.8441652 | 0.7296185 | 0.00e+00 |
| 713 | Mucopolysaccharidoses | 11 | 1.36e-05 | 0.7559480 | -0.2581369 | 0.0000141 | 0.1382052 | 0.7988066 | 0.7170663 | 4.58e-05 |
| 265 | Defective B4GALT7 causes EDS, progeroid type | 18 | 1.00e-07 | 0.3708989 | -0.6386151 | 0.0064379 | 0.0000027 | 0.7385088 | 0.7138342 | 6.00e-07 |
| 617 | L13a-mediated translational silencing of Ceruloplasmin expression | 109 | 0.00e+00 | -0.4660523 | 0.5254805 | 0.0000000 | 0.0000000 | 0.7023778 | 0.7011196 | 0.00e+00 |
| 460 | GTP hydrolysis and joining of the 60S ribosomal subunit | 110 | 0.00e+00 | -0.4418892 | 0.5302721 | 0.0000000 | 0.0000000 | 0.6902568 | 0.6874218 | 0.00e+00 |
| 264 | Defective B3GAT3 causes JDSSDHD | 18 | 9.00e-07 | 0.3764511 | -0.5746291 | 0.0056854 | 0.0000243 | 0.6869600 | 0.6725153 | 4.00e-06 |
| 806 | Other semaphorin interactions | 15 | 8.70e-06 | 0.3594897 | -0.5881920 | 0.0159171 | 0.0000798 | 0.6893495 | 0.6701122 | 3.06e-05 |
| 262 | Defective B3GALT6 causes EDSP2 and SEMDJL1 | 18 | 1.30e-06 | 0.3746746 | -0.5643972 | 0.0059170 | 0.0000338 | 0.6774402 | 0.6640240 | 5.60e-06 |
mitch_report(mres,outfile="mitchmulti_discord.html",overwrite=TRUE)## Note: overwriting existing report## Dataset saved as " /tmp/RtmpMZfPOC/mitchmulti_discord.rds ".##
##
## processing file: mitch.Rmd##
|
| | 0%
|
|.. | 3%
## inline R code fragments
##
##
|
|.... | 6%
## label: checklibraries (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|...... | 9%
## ordinary text without R code
##
##
|
|........ | 12%
## label: peek (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.......... | 15%
## ordinary text without R code
##
##
|
|............ | 18%
## label: metrics (with options)
## List of 1
## $ echo: logi FALSE
##
##
|
|.............. | 21%
## ordinary text without R code
##
##
|
|................ | 24%
## label: scatterplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ message : logi FALSE
## $ warning : logi FALSE##
|
|................... | 26%
## ordinary text without R code
##
##
|
|..................... | 29%
## label: contourplot (with options)
## List of 5
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6.5
## $ warning : logi FALSE
## $ message : logi FALSE##
|
|....................... | 32%
## ordinary text without R code
##
##
|
|......................... | 35%
## label: input_geneset_metrics1 (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|........................... | 38%
## ordinary text without R code
##
##
|
|............................. | 41%
## label: input_geneset_metrics2 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"##
|
|............................... | 44%
## ordinary text without R code
##
##
|
|................................. | 47%
## label: input_geneset_metrics3 (with options)
## List of 5
## $ results : chr "asis"
## $ echo : logi FALSE
## $ message : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7##
|
|................................... | 50%
## ordinary text without R code
##
##
|
|..................................... | 53%
## label: echart1d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|....................................... | 56%
## label: echart2d (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|......................................... | 59%
## ordinary text without R code
##
##
|
|........................................... | 62%
## label: heatmap (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 10
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE##
|
|............................................. | 65%
## ordinary text without R code
##
##
|
|............................................... | 68%
## label: effectsize (with options)
## List of 6
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 7
## $ fig.width : num 7
## $ fig.show : chr "all"
## $ message : logi FALSE
##
##
|
|................................................. | 71%
## ordinary text without R code
##
##
|
|................................................... | 74%
## label: results_table (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|...................................................... | 76%
## ordinary text without R code
##
##
|
|........................................................ | 79%
## label: results_table_complete (with options)
## List of 2
## $ results: chr "asis"
## $ echo : logi FALSE
##
##
|
|.......................................................... | 82%
## ordinary text without R code
##
##
|
|............................................................ | 85%
## label: detailed_geneset_reports1d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 6
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE
##
##
|
|.............................................................. | 88%
## ordinary text without R code
##
##
|
|................................................................ | 91%
## label: detailed_geneset_reports2d (with options)
## List of 7
## $ results : chr "asis"
## $ echo : logi FALSE
## $ fig.height: num 5
## $ fig.width : num 6
## $ out.width : chr "80%"
## $ comment : logi NA
## $ message : logi FALSE##
|
|.................................................................. | 94%
## ordinary text without R code
##
##
|
|.................................................................... | 97%
## label: session_info (with options)
## List of 3
## $ include: logi TRUE
## $ echo : logi TRUE
## $ results: chr "markup"
##
##
|
|......................................................................| 100%
## ordinary text without R code## output file: /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md## /home/mdz/anaconda3/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/mmckenzie/echs1_rnaseq/dge/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpMZfPOC/mitch_report.html --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/lib/R/library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --standalone --section-divs --template /usr/lib/R/library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --include-in-header /tmp/RtmpMZfPOC/rmarkdown-str1793083c26e9.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'##
## Output created: /tmp/RtmpMZfPOC/mitch_report.html## [1] TRUEwrite.table(mres$enrichment_result,file="mitch_multi_discord.tsv",quote=FALSE,sep="\t",row.names=FALSE)
n=40
mx <- as.matrix(mres$enrichment_result[1:n,4:5])
rownames(mx) <- mres$set
rownames(mx) <- mres$enrichment_result$set[1:n]
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2(mx,trace="none",col=colfunc(25),scale="none", margins = c(5,15), cexRow=0.7,
cexCol=0.9, main=paste("Top ranked Reactomes") )
Conclusion
Session information
For reproducibility.
sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] pkgload_1.2.1 GGally_2.1.2
## [3] beeswarm_0.4.0 gtools_3.9.2
## [5] echarts4r_0.4.1 kableExtra_1.3.4
## [7] topconfects_1.8.0 limma_3.48.1
## [9] eulerr_6.1.0 mitch_1.4.0
## [11] MASS_7.3-55 fgsea_1.18.0
## [13] gplots_3.1.1 DESeq2_1.32.0
## [15] SummarizedExperiment_1.22.0 Biobase_2.52.0
## [17] MatrixGenerics_1.4.0 matrixStats_0.59.0
## [19] GenomicRanges_1.44.0 GenomeInfoDb_1.28.1
## [21] IRanges_2.26.0 S4Vectors_0.30.0
## [23] BiocGenerics_0.38.0 reshape2_1.4.4
## [25] forcats_0.5.1 stringr_1.4.0
## [27] dplyr_1.0.7 purrr_0.3.4
## [29] readr_2.0.0 tidyr_1.1.3
## [31] tibble_3.1.3 ggplot2_3.3.5
## [33] tidyverse_1.3.1 zoo_1.8-9
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.2.1 fastmatch_1.1-3
## [4] systemfonts_1.0.2 plyr_1.8.6 polylabelr_0.2.0
## [7] splines_4.1.2 BiocParallel_1.26.1 digest_0.6.27
## [10] htmltools_0.5.1.1 fansi_0.5.0 magrittr_2.0.1
## [13] memoise_2.0.0 tzdb_0.1.2 Biostrings_2.60.1
## [16] annotate_1.70.0 modelr_0.1.8 svglite_2.0.0
## [19] colorspace_2.0-2 blob_1.2.2 rvest_1.0.0
## [22] haven_2.4.1 xfun_0.24 crayon_1.4.1
## [25] RCurl_1.98-1.3 jsonlite_1.7.2 genefilter_1.74.0
## [28] survival_3.2-13 glue_1.4.2 polyclip_1.10-0
## [31] gtable_0.3.0 zlibbioc_1.38.0 XVector_0.32.0
## [34] webshot_0.5.2 DelayedArray_0.18.0 scales_1.1.1
## [37] DBI_1.1.1 Rcpp_1.0.7 viridisLite_0.4.0
## [40] xtable_1.8-4 bit_4.0.4 htmlwidgets_1.5.3
## [43] httr_1.4.2 RColorBrewer_1.1-2 ellipsis_0.3.2
## [46] farver_2.1.0 pkgconfig_2.0.3 reshape_0.8.8
## [49] XML_3.99-0.6 sass_0.4.0 dbplyr_2.1.1
## [52] locfit_1.5-9.4 utf8_1.2.2 labeling_0.4.2
## [55] tidyselect_1.1.1 rlang_0.4.11 later_1.2.0
## [58] AnnotationDbi_1.54.1 munsell_0.5.0 cellranger_1.1.0
## [61] tools_4.1.2 cachem_1.0.5 cli_3.0.1
## [64] generics_0.1.0 RSQLite_2.2.7 broom_0.7.8
## [67] evaluate_0.14 fastmap_1.1.0 yaml_2.2.1
## [70] knitr_1.33 bit64_4.0.5 fs_1.5.0
## [73] caTools_1.18.2 KEGGREST_1.32.0 mime_0.11
## [76] xml2_1.3.2 compiler_4.1.2 rstudioapi_0.13
## [79] png_0.1-7 reprex_2.0.0 geneplotter_1.70.0
## [82] bslib_0.2.5.1 stringi_1.7.3 highr_0.9
## [85] desc_1.3.0 lattice_0.20-45 Matrix_1.4-0
## [88] vctrs_0.3.8 pillar_1.6.1 lifecycle_1.0.0
## [91] jquerylib_0.1.4 data.table_1.14.0 bitops_1.0-7
## [94] httpuv_1.6.1 R6_2.5.0 promises_1.2.0.1
## [97] KernSmooth_2.23-20 gridExtra_2.3 assertthat_0.2.1
## [100] rprojroot_2.0.2 withr_2.4.2 GenomeInfoDbData_1.2.6
## [103] hms_1.1.0 grid_4.1.2 rmarkdown_2.9
## [106] shiny_1.6.0 lubridate_1.7.10