After Kallisto and DESeq2, here I am performing mitch to compare the effect of sedentary and trained rats.
Import the Kallisto transcript counts. We can also include some info out of the Ensembl GTF file including gene name and gene class.
# libraries
library("reshape2")
library("DESeq2")
library("mitch")
library("gplots")
# heres some DESeq2 tables I prepared earlier
skel <- read.table("skeletal_rna.tsv",header=TRUE,row.names=1,sep="\t")
rownames(skel) <- sapply( strsplit(rownames(skel)," ") , "[[", 2)
rownames(skel) <- sapply( strsplit(rownames(skel),"\\.") , "[[", 1)
heart <- read.table("heart_rna.tsv",header=TRUE,row.names=1,sep="\t")
rownames(heart) <- sapply( strsplit(rownames(heart)," ") , "[[", 2)
rownames(heart) <- sapply( strsplit(rownames(heart),"\\.") , "[[", 1)
# genesets and gene table
genesets<-gmt_import("ReactomePathways.gmt")
mart <- read.table("mart_export.txt",header=TRUE,sep="\t")
gt <- mart[,c("Gene.stable.ID","Human.gene.name")]
l <- list("skel"=as.data.frame(skel),"heart"=as.data.frame(heart))
m <- mitch_import(l,"DESeq2",geneTable=gt)## Note: Mean no. genes in input = 32852## Note: no. genes in output = 17470## Note: estimated proportion of input genes in output = 0.532head(res$enrichment_result,25)## set
## 1339 Translation
## 708 Mitochondrial translation elongation
## 707 Mitochondrial translation
## 710 Mitochondrial translation termination
## 709 Mitochondrial translation initiation
## 1299 The citric acid (TCA) cycle and respiratory electron transport
## 662 Macroautophagy
## 686 Metabolism of proteins
## 103 Autophagy
## 1402 Xenobiotics
## 922 Protein localization
## 1066 Respiratory electron transport
## 671 Metabolism
## 272 Defective CFTR causes cystic fibrosis
## 1436 tRNA Aminoacylation
## 380 Extracellular matrix organization
## 1067 Respiratory electron transport, ATP synthesis by chemiosmotic coupling, and heat production by uncoupling proteins.
## 2 ABC transporter disorders
## 1123 Selective autophagy
## 4 ABC-family proteins mediated transport
## 533 Hh mutants that don't undergo autocatalytic processing are degraded by ERAD
## 411 Fatty acids
## 532 Hh mutants abrogate ligand secretion
## 705 Mitochondrial protein import
## 585 Interleukin-1 signaling
## setSize pMANOVA s.skel s.heart p.skel p.heart
## 1339 261 5.957737e-22 0.143023042 0.32118585 7.106377e-05 4.274679e-19
## 708 85 1.131508e-21 0.174218478 0.58527331 5.513769e-03 1.016967e-20
## 707 91 3.596434e-21 0.153504071 0.56295491 1.141470e-02 1.570999e-20
## 710 85 6.428033e-21 0.174429613 0.57356342 5.456992e-03 5.850722e-20
## 709 85 5.598449e-20 0.191220288 0.55276388 2.316823e-03 1.201183e-18
## 1299 164 6.827954e-13 0.081009195 0.32648024 7.371618e-02 5.494891e-13
## 662 108 4.641845e-11 -0.012630820 0.38419793 8.207077e-01 5.260760e-12
## 686 1803 1.013154e-10 0.029691902 0.09181123 3.864540e-02 1.577314e-10
## 103 122 1.766200e-10 -0.044779575 0.34963165 3.933072e-01 2.571702e-11
## 1402 20 3.532813e-10 -0.675710602 -0.49871347 1.670417e-07 1.127604e-04
## 922 156 4.728736e-10 -0.018297035 0.30397514 6.935612e-01 5.749328e-11
## 1066 95 1.466387e-09 0.149201060 0.34353351 1.200514e-02 7.204745e-09
## 671 1930 8.450875e-09 -0.026260061 0.08149145 5.949212e-02 4.893389e-09
## 272 58 1.484545e-08 0.059867551 0.44986216 4.304944e-01 3.089215e-09
## 1436 41 1.987553e-08 0.211073778 0.48785666 1.937746e-02 6.440694e-08
## 380 271 2.111213e-08 0.125158096 -0.17227016 3.982399e-04 1.087928e-06
## 1067 116 2.807109e-08 0.123549972 0.28826803 2.161081e-02 8.252947e-08
## 2 74 2.906455e-08 0.032411924 0.39363818 6.298982e-01 4.767368e-09
## 1123 56 4.621963e-08 -0.112999188 0.43759844 1.436788e-01 1.470053e-08
## 4 97 4.800773e-08 0.072211234 0.33116680 2.193173e-01 1.743128e-08
## 533 53 1.655431e-07 0.055442471 0.43841140 4.851793e-01 3.355831e-08
## 411 14 1.933391e-07 -0.708732978 -0.46281671 4.383211e-06 2.714328e-03
## 532 56 2.124180e-07 0.045406816 0.42430454 5.568300e-01 3.956725e-08
## 705 62 2.593500e-07 0.129682682 0.37910082 7.748860e-02 2.434610e-07
## 585 94 3.536827e-07 0.005196089 0.32501849 9.306589e-01 5.174549e-08
## s.dist SD p.adjustMANOVA
## 1339 0.35159059 0.12598013 8.158175e-19
## 708 0.61065287 0.29065966 8.158175e-19
## 707 0.58350813 0.28952547 1.728686e-18
## 710 0.59950036 0.28223022 2.317306e-18
## 709 0.58490435 0.25564992 1.614593e-17
## 1299 0.33638050 0.17357424 1.640985e-10
## 662 0.38440550 0.28060030 9.562201e-09
## 686 0.09649306 0.04392500 1.826211e-08
## 103 0.35248759 0.27889085 2.829844e-08
## 1402 0.83982137 0.12515587 5.094316e-08
## 922 0.30452532 0.22788084 6.198943e-08
## 1066 0.37453468 0.13741379 1.762109e-07
## 671 0.08561804 0.07619183 9.373971e-07
## 272 0.45382826 0.27576784 1.529082e-06
## 1436 0.53156021 0.19571505 1.902731e-06
## 380 0.21293557 0.21031354 1.902731e-06
## 1067 0.31362885 0.11647326 2.328393e-06
## 2 0.39497031 0.25542553 2.328393e-06
## 1123 0.45195267 0.38933132 3.461358e-06
## 4 0.33894824 0.18310924 3.461358e-06
## 533 0.44190318 0.27079993 1.136729e-05
## 411 0.84646426 0.17388906 1.267250e-05
## 532 0.42672722 0.26792115 1.331769e-05
## 705 0.40066823 0.17636526 1.558261e-05
## 585 0.32506002 0.22614859 2.040042e-05unlink("2d_mitch1.html")
mitch_report(res,outfile="2d_mitch1.html")## Dataset saved as " /tmp/RtmptTLcro/2d_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##
|
|....................... | 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/adam_trewin/integrate/mitch.knit.md## /usr/bin/pandoc +RTS -K512m -RTS /mnt/bfx6/bfx/adam_trewin/integrate/mitch.utf8.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmptTLcro/mitch_report.html --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/RtmptTLcro/rmarkdown-str4e2eea702cf.html --mathjax --variable 'mathjax-url:https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'##
## Output created: /tmp/RtmptTLcro/mitch_report.html## [1] TRUEsessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## 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] ggplot2_3.3.2 beeswarm_0.2.3
## [5] gtools_3.8.2 tibble_3.0.4
## [7] dplyr_1.0.2 echarts4r_0.3.3
## [9] gplots_3.1.0 mitch_1.2.2
## [11] DESeq2_1.30.0 SummarizedExperiment_1.20.0
## [13] Biobase_2.50.0 MatrixGenerics_1.2.0
## [15] matrixStats_0.57.0 GenomicRanges_1.42.0
## [17] GenomeInfoDb_1.26.0 IRanges_2.24.0
## [19] S4Vectors_0.28.0 BiocGenerics_0.36.0
## [21] reshape2_1.4.4
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-6 bit64_4.0.5 RColorBrewer_1.1-2
## [4] httr_1.4.2 rprojroot_1.3-2 tools_4.0.3
## [7] backports_1.2.0 R6_2.5.0 KernSmooth_2.23-18
## [10] DBI_1.1.0 colorspace_2.0-0 withr_2.3.0
## [13] tidyselect_1.1.0 gridExtra_2.3 bit_4.0.4
## [16] compiler_4.0.3 desc_1.2.0 DelayedArray_0.16.0
## [19] labeling_0.4.2 caTools_1.18.0 scales_1.1.1
## [22] genefilter_1.72.0 stringr_1.4.0 digest_0.6.27
## [25] rmarkdown_2.5 XVector_0.30.0 pkgconfig_2.0.3
## [28] htmltools_0.5.0 fastmap_1.0.1 highr_0.8
## [31] htmlwidgets_1.5.2 rlang_0.4.8 RSQLite_2.2.1
## [34] shiny_1.5.0 farver_2.0.3 generics_0.1.0
## [37] jsonlite_1.7.1 BiocParallel_1.24.1 RCurl_1.98-1.2
## [40] magrittr_1.5 GenomeInfoDbData_1.2.4 Matrix_1.2-18
## [43] Rcpp_1.0.5 munsell_0.5.0 lifecycle_0.2.0
## [46] stringi_1.5.3 yaml_2.2.1 MASS_7.3-53
## [49] zlibbioc_1.36.0 plyr_1.8.6 grid_4.0.3
## [52] blob_1.2.1 promises_1.1.1 crayon_1.3.4
## [55] lattice_0.20-41 splines_4.0.3 annotate_1.68.0
## [58] locfit_1.5-9.4 knitr_1.30 pillar_1.4.6
## [61] geneplotter_1.68.0 XML_3.99-0.5 glue_1.4.2
## [64] evaluate_0.14 vctrs_0.3.4 httpuv_1.5.4
## [67] testthat_3.0.0 gtable_0.3.0 purrr_0.3.4
## [70] reshape_0.8.8 assertthat_0.2.1 xfun_0.19
## [73] mime_0.9 xtable_1.8-4 later_1.1.0.1
## [76] survival_3.2-7 AnnotationDbi_1.52.0 memoise_1.1.0
## [79] ellipsis_0.3.1