Analysis of mir-101 RNA-seq data
Mark Ziemann
2021-03-02
introduction
run deseq2 with Matt’s counts from AGRF
library("tidyverse")
library("DESeq2")
library("fgsea")
library("gplots")
library("mitch")
library("vioplot")
library("kableExtra")import data
read counts from AGRF and genrate MDS plot.
MM1 samples are the treatment group. MMN are the ctrl group. Matt please check this.
# counts
counts <- read.table("CAGRF13809_cnt.tsv",header=T,row.names=1)
head(counts)## MM1_1 MM1_2 MM1_3 MM1_4 MMN_1 MMN_2 MMN_3 MMN_4
## CCDC80 1592 1779 1706 1583 420 431 286 342
## STK25 2159 2576 2494 2196 560 692 459 762
## TLR4 480 508 558 495 1523 1447 1480 1365
## LUM 39 38 35 43 1 4 1 0
## SH2B3 234 228 320 155 765 1005 790 995
## TFB2M 325 340 375 418 867 909 781 1023colnames(counts) <- gsub("MMN","CTRL",colnames(counts))
colnames(counts) <- gsub("MM1","M101",colnames(counts))
head(counts)## M101_1 M101_2 M101_3 M101_4 CTRL_1 CTRL_2 CTRL_3 CTRL_4
## CCDC80 1592 1779 1706 1583 420 431 286 342
## STK25 2159 2576 2494 2196 560 692 459 762
## TLR4 480 508 558 495 1523 1447 1480 1365
## LUM 39 38 35 43 1 4 1 0
## SH2B3 234 228 320 155 765 1005 790 995
## TFB2M 325 340 375 418 867 909 781 1023# sampleshet
des<-as.data.frame(colnames(counts))
des$grp<-as.numeric(grepl("M101",colnames(counts)))
rownames(des)<-des[,1]
des[,1]=NULL
colours = c('lightblue', 'pink')
mds <- cmdscale(dist(t(counts)))
XMAX=max(mds[,1])*1.1
XMIN=min(mds[,1])*1.1
plot( mds , xlab="Coordinate 1", ylab="Coordinate 2", pch=19, cex=1.5,
col= colours[as.factor(des$grp)], type = "p" , xlim=c(XMIN,XMAX),main="MDS plot",bty="n" )
legend('topright', col=colours, legend=c("Ctrl","miR-101"), pch = 16, cex = 1)
text(mds*1.05, labels=colnames(counts) )
pdf("mds.pdf")
plot( mds , xlab="Coordinate 1", ylab="Coordinate 2", pch=19, cex=1.5,
col= colours[as.factor(des$grp)], type = "p" , xlim=c(XMIN,XMAX),main="MDS plot",bty="n" )
legend('topright', col=colours, legend=c("Ctrl","miR-101"), pch = 16, cex = 1)
text(mds*1.05, labels=colnames(counts) )
dev.off()## X11cairo
## 2find stable genes
There is a problem with the qRT-PCR validation. So here I will identify some stable genes based on low coefficient of variation. We can also look at expression of housekeeping genes ACTB and B2M. ACTB and B2M are actually not very stable genes according to this data and there thousands with more stable expression.
head(counts)## M101_1 M101_2 M101_3 M101_4 CTRL_1 CTRL_2 CTRL_3 CTRL_4
## CCDC80 1592 1779 1706 1583 420 431 286 342
## STK25 2159 2576 2494 2196 560 692 459 762
## TLR4 480 508 558 495 1523 1447 1480 1365
## LUM 39 38 35 43 1 4 1 0
## SH2B3 234 228 320 155 765 1005 790 995
## TFB2M 325 340 375 418 867 909 781 1023x <- counts/colSums(counts) * 1e6
xsd <- apply(x,1,sd)
xmean <- rowMeans(x)
xcv <- apply(x,1,sd) / rowMeans(x)
x <- data.frame(x,xmean,xsd,xcv)
head( x[order(x$xcv),] , 20) %>% kbl() %>% kable_styling()| M101_1 | M101_2 | M101_3 | M101_4 | CTRL_1 | CTRL_2 | CTRL_3 | CTRL_4 | xmean | xsd | xcv | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| TFG | 67.725154 | 70.504046 | 70.758367 | 72.313140 | 70.921621 | 67.073566 | 69.791166 | 74.843831 | 70.491361 | 2.4594700 | 0.0348904 |
| EFR3A | 160.872212 | 169.834365 | 168.649258 | 165.033339 | 162.570335 | 172.240821 | 178.168550 | 174.456124 | 168.978126 | 5.9556190 | 0.0352449 |
| FBXL3 | 81.609809 | 88.219660 | 85.520904 | 88.520329 | 86.903958 | 92.571761 | 85.673620 | 89.866441 | 87.360810 | 3.2738990 | 0.0374756 |
| AHCTF1 | 95.285425 | 101.403357 | 106.827454 | 95.594354 | 98.818307 | 102.574179 | 101.658383 | 100.588834 | 100.343787 | 3.7870968 | 0.0377412 |
| DPH6 | 17.523046 | 17.894622 | 19.111919 | 17.212074 | 16.911130 | 17.383347 | 18.230678 | 17.162471 | 17.678661 | 0.7161681 | 0.0405103 |
| UBL3 | 48.796073 | 50.484379 | 48.054603 | 49.752302 | 46.848226 | 48.743538 | 53.145133 | 47.060078 | 49.110541 | 2.0440952 | 0.0416223 |
| LRRC8B | 31.130327 | 32.273960 | 33.491268 | 34.062356 | 34.509606 | 34.921036 | 34.783536 | 32.114509 | 33.410825 | 1.4130476 | 0.0422931 |
| COX16 | 4.060896 | 4.033391 | 3.965585 | 4.017804 | 4.450297 | 3.976583 | 3.965585 | 4.265816 | 4.091995 | 0.1747566 | 0.0427069 |
| BRCA2 | 207.634425 | 217.556429 | 202.378272 | 229.134684 | 207.918467 | 209.845569 | 201.336619 | 207.383855 | 210.398540 | 9.0388586 | 0.0429607 |
| C11orf58 | 169.222556 | 170.197740 | 179.387668 | 159.174479 | 161.879565 | 167.754982 | 164.021024 | 156.595148 | 166.029145 | 7.2332764 | 0.0435663 |
| DNA2 | 58.243266 | 62.375540 | 60.475180 | 59.523021 | 59.578356 | 65.102340 | 60.199792 | 56.298857 | 60.224544 | 2.6384323 | 0.0438099 |
| FAM199X | 187.293973 | 189.520431 | 205.901301 | 186.144283 | 204.855983 | 200.210544 | 207.462791 | 196.982305 | 197.296451 | 8.6851432 | 0.0440208 |
| RBM27 | 60.775819 | 58.541095 | 63.159578 | 65.177969 | 59.086180 | 57.166652 | 63.105734 | 61.731777 | 61.093100 | 2.7116998 | 0.0443864 |
| ERCC6L2 | 86.097328 | 85.655338 | 88.885275 | 86.182673 | 87.497284 | 82.508815 | 79.669060 | 92.342214 | 86.104748 | 3.8384033 | 0.0445783 |
| MYL12B | 398.449171 | 411.545181 | 429.936438 | 428.469909 | 436.894132 | 411.495236 | 442.531932 | 390.138218 | 418.682527 | 18.7070084 | 0.0446807 |
| GHITM | 356.719697 | 374.236413 | 380.015272 | 359.085985 | 406.256621 | 382.776974 | 358.049959 | 363.362030 | 372.562869 | 17.0119639 | 0.0456620 |
| ATAD1 | 180.842622 | 180.357477 | 191.470980 | 189.390695 | 178.333764 | 174.605178 | 186.463443 | 165.966582 | 180.928843 | 8.3177754 | 0.0459726 |
| FGGY | 7.328274 | 7.655783 | 6.894750 | 6.953590 | 6.703382 | 7.105007 | 6.944352 | 6.619817 | 7.025620 | 0.3364606 | 0.0478905 |
| METTL23 | 34.099996 | 36.499100 | 39.252628 | 36.359817 | 38.349695 | 35.379183 | 37.206538 | 34.811528 | 36.494811 | 1.7478211 | 0.0478923 |
| AK3 | 120.503952 | 113.424649 | 129.897433 | 114.435114 | 123.734621 | 121.365872 | 116.585122 | 114.231493 | 119.272282 | 5.7187639 | 0.0479471 |
x[which(rownames(x) == "ACTB"),]## M101_1 M101_2 M101_3 M101_4 CTRL_1 CTRL_2 CTRL_3 CTRL_4
## ACTB 288.4449 278.3424 299.6806 389.8328 329.7436 273.8474 304.8007 688.0869
## xmean xsd xcv
## ACTB 356.5974 139.0014 0.3897991x[which(rownames(x) == "B2M"),]## M101_1 M101_2 M101_3 M101_4 CTRL_1 CTRL_2 CTRL_3 CTRL_4
## B2M 195.8183 183.7388 140.9208 151.032 188.1492 171.7319 129.4626 84.88942
## xmean xsd xcv
## B2M 155.7179 37.12599 0.2384184hist(xcv,breaks=50,xlim=c(0,1),main="histogram of CV values")
DESeq2
Now run DESeq2 to identify differential expression. Then make some charts.
dds <- DESeqDataSetFromMatrix(countData = counts , colData = des, design = ~ grp)## the design formula contains one or more numeric variables with integer values,
## specifying a model with increasing fold change for higher values.
## did you mean for this to be a factor? if so, first convert
## this variable to a factor using the factor() functionres <- 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)
#stick on the normalised expression values to the table
zz<-cbind(as.data.frame(z),assay(vsd))
#sort by p-value
mm1<-as.data.frame(zz[order(zz$pvalue),])
#some plots
sig<-subset(zz,padj<0.05)
SIG=nrow(sig)
DN=nrow(subset(sig,log2FoldChange<0))
UP=nrow(subset(sig,log2FoldChange>0))
HEADER=paste("Ctrl vs miR-101:", SIG , "DGEs,", UP ,"upregulated,", DN, "downregulated")
plot(log2(zz$baseMean),zz$log2FoldChange,cex=0.6, xlab="log2 base mean",
ylab="log2 fold change" ,pch=19,col="#838383")
points(log2(sig$baseMean),sig$log2FoldChange,cex=0.6,pch=19,col="red")
mtext(HEADER)
top<-head(sig,20)
text(log2(top$baseMean)+1, top$log2FoldChange, labels = rownames(top),cex=0.7)
#volcano plot
plot(zz$log2FoldChange, -log2(zz$pvalue) ,cex=0.6, xlim=c(-4,6),
xlab="log2 fold change", ylab="-log2 p-value" ,pch=19,col="#838383")
points(sig$log2FoldChange, -log2(sig$pvalue),cex=0.6,pch=19,col="red")
text(top$log2FoldChange+0.5, -log2(top$pvalue), labels = rownames(top),cex=0.7)
mtext(HEADER)
pdf("volcanoplot.pdf")
#volcano plot
plot(zz$log2FoldChange, -log2(zz$pvalue) ,cex=0.6, xlim=c(-4,6),
xlab="log2 fold change", ylab="-log2 p-value" ,pch=19,col="#838383")
points(sig$log2FoldChange, -log2(sig$pvalue),cex=0.6,pch=19,col="red")
text(top$log2FoldChange+0.5, -log2(top$pvalue), labels = rownames(top),cex=0.8)
mtext(HEADER)
dev.off()## X11cairo
## 2# top N genes
colfunc <- colorRampPalette(c("blue", "white", "red"))
heatmap.2( as.matrix(zz[1:50,c(7:ncol(zz))]), col=colfunc(25),scale="row",
trace="none",margins = c(6,10), cexRow=0.7, main="Top 50 genes")
pdf("heatmap.pdf")
heatmap.2( as.matrix(zz[1:50,c(7:ncol(zz))]), col=colfunc(25),scale="row",
trace="none",margins = c(6,10), cexRow=0.5, main="Top 50 genes")
dev.off()## X11cairo
## 2#output DGE table
write.table(zz,file="CAGRF13809_deseq.tsv",quote=F,sep="\t")Run the mitch package with gene sets from Reactome. Check out the mitchreport.html file for the full report.
# gene set downloaded 2021-Jan-08
#download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip",
# destfile="ReactomePathways.gmt.zip")
#unzip("ReactomePathways.gmt.zip")
genesets <- gmt_import("ReactomePathways.gmt")
y <- mitch_import(mm1, DEtype="deseq2")## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 15982## Note: no. genes in output = 15982## Note: estimated proportion of input genes in output = 1res <- mitch_calc(y, genesets, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(res$enrichment_result,30)## set
## 522 Interleukin-35 Signalling
## 755 Peptide chain elongation
## 1253 Viral mRNA Translation
## 1000 Selenocysteine synthesis
## 320 Eukaryotic Translation Elongation
## 322 Eukaryotic Translation Termination
## 999 Selenoamino acid metabolism
## 948 Response of EIF2AK4 (GCN2) to amino acid deficiency
## 695 Nonsense Mediated Decay (NMD) independent of the Exon Junction Complex (EJC)
## 526 Interleukin-6 signaling
## 360 Formation of a pool of free 40S subunits
## 31 Activation of BAD and translocation to mitochondria
## 1050 Signaling by Leptin
## 156 Chemokine receptors bind chemokines
## 358 Formation of Senescence-Associated Heterochromatin Foci (SAHF)
## 546 L13a-mediated translational silencing of Ceruloplasmin expression
## 834 RHO GTPases activate KTN1
## 245 Degradation of cysteine and homocysteine
## 400 GTP hydrolysis and joining of the 60S ribosomal subunit
## 1115 Synthesis of PIPs at the Golgi membrane
## 833 RHO GTPases activate IQGAPs
## 572 MAPK3 (ERK1) activation
## 1116 Synthesis of PIPs at the early endosome membrane
## 1298 tRNA processing in the mitochondrion
## 980 SRP-dependent cotranslational protein targeting to membrane
## 127 Cap-dependent Translation Initiation
## 321 Eukaryotic Translation Initiation
## 306 Early Phase of HIV Life Cycle
## 491 Inhibition of replication initiation of damaged DNA by RB1/E2F1
## 527 Interleukin-7 signaling
## setSize pANOVA s.dist p.adjustANOVA
## 522 10 2.520035e-04 0.6683822 3.276045e-03
## 755 85 1.380194e-24 0.6416579 1.794253e-21
## 1253 85 4.653191e-24 0.6342647 1.928667e-21
## 1000 88 3.640477e-24 0.6248928 1.928667e-21
## 320 89 7.417951e-24 0.6171306 1.928667e-21
## 322 89 2.092458e-22 0.5967343 3.885994e-20
## 999 101 6.140168e-24 0.5805948 1.928667e-21
## 948 97 9.129945e-23 0.5766557 1.978155e-20
## 695 91 1.493911e-20 0.5633550 1.787329e-18
## 526 11 1.284997e-03 0.5605558 1.084737e-02
## 360 97 1.393878e-21 0.5603438 2.265051e-19
## 31 12 8.120316e-04 -0.5582655 7.705409e-03
## 1050 10 2.687708e-03 0.5481092 1.792302e-02
## 156 13 9.708503e-04 0.5284036 8.888066e-03
## 358 11 2.568127e-03 -0.5250255 1.765531e-02
## 546 107 8.216220e-21 0.5233286 1.186787e-18
## 834 10 4.533943e-03 -0.5183571 2.529668e-02
## 245 11 2.964013e-03 -0.5174094 1.870494e-02
## 400 108 1.512355e-20 0.5173239 1.787329e-18
## 1115 15 5.864811e-04 -0.5126615 6.395300e-03
## 833 10 5.312492e-03 -0.5090533 2.796048e-02
## 572 10 5.482446e-03 0.5071876 2.861483e-02
## 1116 15 8.035164e-04 -0.4998100 7.689560e-03
## 1298 24 3.005470e-05 -0.4920729 5.502973e-04
## 980 108 2.112169e-18 0.4872795 2.288183e-16
## 127 115 2.624925e-18 0.4709983 2.437431e-16
## 321 115 2.624925e-18 0.4709983 2.437431e-16
## 306 12 4.768027e-03 -0.4705385 2.615083e-02
## 491 12 4.979772e-03 -0.4682112 2.682830e-02
## 527 17 8.953125e-04 0.4653083 8.434104e-03write.table(res$enrichment_result,file="CAGRF13809_mitch.tsv",quote=F,sep="\t")
unlink("mitchreport.html")
mitch_report(res, "mitchreport.html")## Dataset saved as " /tmp/Rtmp98ceqO/mitchreport.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/deseq2/mitch.knit.md## /usr/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/mmckenzie/deseq2/mitch.utf8.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/Rtmp98ceqO/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --email-obfuscation none --self-contained --standalone --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable 'theme:bootstrap' --include-in-header /tmp/Rtmp98ceqO/rmarkdown-str1d301076e5f1b5.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'##
## Output created: /tmp/Rtmp98ceqO/mitch_report.html## [1] TRUEtop <- res$enrichment_result
top <- head(subset(top,p.adjustANOVA<0.05),25)
top <- top[order(-top$s.dist),]
par(mar=c(5, 29, 4, 2))
barplot(top$s.dist,horiz=TRUE,las=1,names.arg=top$set,cex.names=0.8, xlab="enrichment score")
pdf("reactome_barplot.pdf")
par(mar=c(5, 29, 4, 2))
barplot(top$s.dist,horiz=TRUE,las=1,names.arg=top$set,cex.names=0.8, xlab="enrichment score")
dev.off()## X11cairo
## 2par(mar=c(5, 5, 4, 2))
pdf("reactome_volplot.pdf")
plot(res$enrichment_result$s.dist,-log10(res$enrichment_result$pANOVA),
col="gray",pch=19,
xlab="enrichment score", ylab="-log10 p-value")
grid()
sig <- subset(res$enrichment_result,p.adjustANOVA<0.05)
points(sig$s.dist,-log10(sig$pANOVA), col="red",pch=19)
top <- head(sig,5)
mtext("Reactome gene set enrichment")
dev.off()## X11cairo
## 2plot(res$enrichment_result$s.dist,-log10(res$enrichment_result$pANOVA),
col="gray",pch=19,
xlab="enrichment score", ylab="-log10 p-value")
grid()
sig <- subset(res$enrichment_result,p.adjustANOVA<0.05)
points(sig$s.dist,-log10(sig$pANOVA), col="red",pch=19)
top <- head(sig,5)
mtext("Reactome gene set enrichment")
Other genesets to look at.
tRNA processing in the mitochondrion
Complex I biogenesis
The citric acid (TCA) cycle and respiratory electron transport, and/or
Respiratory electron transport, ATP synthesis by chemiosmotic coupling, and heat production by uncoupling proteins.
Regulation of lipid metabolism by PPARalpha
mysets <- c("tRNA processing in the mitochondrion",
"Complex I biogenesis",
"The citric acid (TCA) cycle and respiratory electron transport",
"Respiratory electron transport, ATP synthesis by chemiosmotic coupling, and heat production by uncoupling proteins.",
"Regulation of lipid metabolism by PPARalpha")
mysets <- genesets[which(names(genesets) %in% mysets)]
res <- mitch_calc(y, mysets, priority="effect")## Note: Enrichments with large effect sizes may not be
## statistically significant.head(res$enrichment_result)## set
## 5 tRNA processing in the mitochondrion
## 1 Complex I biogenesis
## 4 The citric acid (TCA) cycle and respiratory electron transport
## 3 Respiratory electron transport, ATP synthesis by chemiosmotic coupling, and heat production by uncoupling proteins.
## 2 Regulation of lipid metabolism by PPARalpha
## setSize pANOVA s.dist p.adjustANOVA
## 5 24 3.005470e-05 -0.4920729 5.009116e-05
## 1 52 7.459594e-05 -0.3176083 9.324493e-05
## 4 147 1.557203e-08 -0.2704252 7.786013e-08
## 3 103 4.323334e-06 -0.2621922 1.080834e-05
## 2 106 5.368419e-04 -0.1947463 5.368419e-04unlink("mygeneset_report.html")
mitch_report(res, "mygeneset_report.html")## Dataset saved as " /tmp/Rtmp98ceqO/mygeneset_report.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/deseq2/mitch.knit.md## /usr/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/mmckenzie/deseq2/mitch.utf8.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/Rtmp98ceqO/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --email-obfuscation none --self-contained --standalone --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable 'theme:bootstrap' --include-in-header /tmp/Rtmp98ceqO/rmarkdown-str1d30101cb8e2f.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'##
## Output created: /tmp/Rtmp98ceqO/mitch_report.html## [1] TRUEvioplot(res$detailed_sets,horizontal=TRUE,las=2,side="right",ylim=c(-8146,7835))
grid()
pdf("reactome_vioplot.pdf")
par(mar=c(5, 29, 4, 2))
vioplot(res$detailed_sets,horizontal=TRUE,las=2,side="right",ylim=c(-8146,7835))
grid()
dev.off()## X11cairo
## 2Starbase
starbase <- gmt_import("starBaseV3.gmt")
res <- mitch_calc(y, starbase, priority="significance",resrows=10)## Note: When prioritising by significance (ie: small
## p-values), large effect sizes might be missed.head(res$enrichment_result,10)## set setSize pANOVA s.dist p.adjustANOVA
## 10 hsa-miR-101-3p 849 3.237607e-100 -0.4297220 1.596140e-97
## 73 hsa-miR-144-3p 888 1.856943e-40 -0.2647680 4.577364e-38
## 294 hsa-miR-423-5p 342 1.283180e-24 -0.3227285 2.108693e-22
## 451 hsa-miR-654-5p 212 2.405229e-21 -0.3781535 2.964445e-19
## 53 hsa-miR-1321 344 3.186886e-21 -0.2971927 3.142270e-19
## 345 hsa-miR-491-5p 377 5.191093e-21 -0.2826543 4.265348e-19
## 338 hsa-miR-485-5p 370 2.224621e-20 -0.2805853 1.454175e-18
## 100 hsa-miR-185-5p 509 2.359716e-20 -0.2401340 1.454175e-18
## 409 hsa-miR-541-3p 237 5.051484e-20 -0.3457855 2.767091e-18
## 351 hsa-miR-497-5p 1139 1.284353e-19 -0.1606638 5.489059e-18sset <- subset(res$enrichment_result,setSize>49)
head( sset[order(sset$s.dist),] )## set setSize pANOVA s.dist p.adjustANOVA
## 10 hsa-miR-101-3p 849 3.237607e-100 -0.4297220 1.596140e-97
## 451 hsa-miR-654-5p 212 2.405229e-21 -0.3781535 2.964445e-19
## 468 hsa-miR-766-5p 80 9.030610e-09 -0.3718007 8.400172e-08
## 207 hsa-miR-3184-5p 94 7.596193e-10 -0.3673079 8.141137e-09
## 380 hsa-miR-5194 54 4.235411e-06 -0.3618901 1.711523e-05
## 223 hsa-miR-331-3p 79 3.712900e-08 -0.3582677 2.732029e-07unlink("starbase_report.html")
mitch_report(res, "starbase_report.html")## Dataset saved as " /tmp/Rtmp98ceqO/starbase_report.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/deseq2/mitch.knit.md## /usr/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/mmckenzie/deseq2/mitch.utf8.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/Rtmp98ceqO/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --email-obfuscation none --self-contained --standalone --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable 'theme:bootstrap' --include-in-header /tmp/Rtmp98ceqO/rmarkdown-str1d301046c80429.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'##
## Output created: /tmp/Rtmp98ceqO/mitch_report.html## [1] TRUEpar(mar=c(5, 20, 4, 2))
vioplot(rev(res$detailed),horizontal=TRUE,las=2,side="right",ylim=c(-8146,7835))
grid()
pdf("starbase_vioplot.pdf")
par(mar=c(5, 20, 4, 2))
vioplot(rev(res$detailed),horizontal=TRUE,las=2,side="right",ylim=c(-8146,7835))
grid()
dev.off()## X11cairo
## 2par(mar=c(5, 5, 4, 2))
pdf("starbase_volplot.pdf")
plot(res$enrichment_result$s.dist,-log10(res$enrichment_result$pANOVA),
col="gray",pch=19,
xlab="enrichment score", ylab="-log10 p-value")
grid()
sig <- subset(res$enrichment_result,p.adjustANOVA<0.05)
points(sig$s.dist,-log10(sig$pANOVA), col="red",pch=19)
top <- head(sig,5)
text(top$s.dist,-log10(top$pANOVA),labels=top$set)
mtext("starBase miR target gene set enrichment")
dev.off()## X11cairo
## 2plot(res$enrichment_result$s.dist,-log10(res$enrichment_result$pANOVA),
col="gray",pch=19,
xlab="enrichment score", ylab="-log10 p-value")
grid()
sig <- subset(res$enrichment_result,p.adjustANOVA<0.05)
points(sig$s.dist,-log10(sig$pANOVA), col="red",pch=19)
top <- head(sig,5)
text(top$s.dist,-log10(top$pANOVA),labels=top$set)
mtext("starBase miR target gene set enrichment")
session information
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 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.1.0 GGally_2.0.0
## [3] reshape2_1.4.4 beeswarm_0.2.3
## [5] gtools_3.8.2 echarts4r_0.3.3
## [7] kableExtra_1.3.4 vioplot_0.3.5
## [9] zoo_1.8-8 sm_2.2-5.6
## [11] mitch_1.2.2 gplots_3.1.0
## [13] fgsea_1.16.0 DESeq2_1.30.0
## [15] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [17] MatrixGenerics_1.2.0 matrixStats_0.57.0
## [19] GenomicRanges_1.42.0 GenomeInfoDb_1.26.0
## [21] IRanges_2.24.0 S4Vectors_0.28.0
## [23] BiocGenerics_0.36.0 forcats_0.5.0
## [25] stringr_1.4.0 dplyr_1.0.2
## [27] purrr_0.3.4 readr_1.4.0
## [29] tidyr_1.1.2 tibble_3.0.4
## [31] ggplot2_3.3.2 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-0 ellipsis_0.3.1 rprojroot_1.3-2
## [4] XVector_0.30.0 fs_1.5.0 rstudioapi_0.12
## [7] bit64_4.0.5 AnnotationDbi_1.52.0 fansi_0.4.1
## [10] lubridate_1.7.9 xml2_1.3.2 splines_4.0.3
## [13] geneplotter_1.68.0 knitr_1.30 jsonlite_1.7.1
## [16] broom_0.7.2 annotate_1.68.0 dbplyr_2.0.0
## [19] shiny_1.5.0 compiler_4.0.3 httr_1.4.2
## [22] backports_1.2.0 assertthat_0.2.1 Matrix_1.2-18
## [25] fastmap_1.0.1 cli_2.1.0 later_1.1.0.1
## [28] htmltools_0.5.0 tools_4.0.3 gtable_0.3.0
## [31] glue_1.4.2 GenomeInfoDbData_1.2.4 fastmatch_1.1-0
## [34] Rcpp_1.0.5 cellranger_1.1.0 vctrs_0.3.4
## [37] svglite_2.0.0 xfun_0.19 ps_1.4.0
## [40] testthat_3.0.0 rvest_0.3.6 mime_0.9
## [43] lifecycle_0.2.0 XML_3.99-0.5 zlibbioc_1.36.0
## [46] MASS_7.3-53 scales_1.1.1 hms_0.5.3
## [49] promises_1.1.1 RColorBrewer_1.1-2 yaml_2.2.1
## [52] memoise_1.1.0 gridExtra_2.3 reshape_0.8.8
## [55] stringi_1.5.3 RSQLite_2.2.1 highr_0.8
## [58] genefilter_1.72.0 desc_1.2.0 caTools_1.18.0
## [61] BiocParallel_1.24.1 systemfonts_1.0.1 rlang_0.4.8
## [64] pkgconfig_2.0.3 bitops_1.0-6 evaluate_0.14
## [67] lattice_0.20-41 htmlwidgets_1.5.2 bit_4.0.4
## [70] tidyselect_1.1.0 plyr_1.8.6 magrittr_1.5
## [73] R6_2.5.0 generics_0.1.0 DelayedArray_0.16.0
## [76] DBI_1.1.0 pillar_1.4.6 haven_2.3.1
## [79] withr_2.3.0 survival_3.2-7 RCurl_1.98-1.2
## [82] modelr_0.1.8 crayon_1.3.4 KernSmooth_2.23-18
## [85] rmarkdown_2.5 locfit_1.5-9.4 grid_4.0.3
## [88] readxl_1.3.1 data.table_1.13.2 blob_1.2.1
## [91] webshot_0.5.2 reprex_0.3.0 digest_0.6.27
## [94] xtable_1.8-4 httpuv_1.5.4 munsell_0.5.0
## [97] viridisLite_0.3.0 tcltk_4.0.3