## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", crop = NULL ) ## ----getPackage, eval=FALSE--------------------------------------------------- # ## install MIRit from Bioconductor # if (!requireNamespace("BiocManager", quietly = TRUE)) { # install.packages("BiocManager") # } # BiocManager::install("MIRit") ## ----getPackageDevel, eval=FALSE---------------------------------------------- # ## install the development version from GitHub # if (!requireNamespace("BiocManager", quietly = TRUE)) { # install.packages("BiocManager") # } # BiocManager::install("jacopo-ronchi/MIRit") ## ----loadPackage, message=FALSE----------------------------------------------- ## load MIRit library(MIRit) ## ----example------------------------------------------------------------------ ## load count matrix for genes data(geneCounts, package = "MIRit") ## load count matrix for miRNAs data(mirnaCounts, package = "MIRit") ## ----geneExpr, echo=FALSE----------------------------------------------------- ## print a table for gene expression matrix knitr::kable(geneCounts[seq(5), seq(5)], digits = 2, caption = "Gene expression matrix. An expression matrix containing normalized values, with row names according to hgnc.") ## ----mirnaExpr, echo=FALSE---------------------------------------------------- ## print a table for miRNA expression matrix knitr::kable(mirnaCounts[seq(5), seq(5)], digits = 2, caption = "MiRNA expression matrix. An expression matrix containing normalized values, with row names according to miRBase.") ## ----colnames----------------------------------------------------------------- ## print sample names in geneCounts colnames(geneCounts) ## print sample names in mirnaCounts colnames(mirnaCounts) ## check if samples in geneCounts are equal to those in mirnaCounts identical(colnames(geneCounts), colnames(mirnaCounts)) ## ----metadata----------------------------------------------------------------- ## create a data.frame containing sample metadata meta <- data.frame(primary = colnames(mirnaCounts), mirnaCol = colnames(mirnaCounts), geneCol = colnames(geneCounts), disease = c(rep("PTC", 8), rep("NTH", 8)), patient = c(rep(paste("Sample_", seq(8), sep = ""), 2))) ## ----mirnaObj----------------------------------------------------------------- ## create the MirnaExperiment object experiment <- MirnaExperiment(mirnaExpr = mirnaCounts, geneExpr = geneCounts, samplesMetadata = meta, pairedSamples = TRUE) ## ----mds, fig.wide=TRUE, fig.cap="MDS plots for miRNAs and genes. Both plots show that the main source of variability is given by the disease condition."---- geneMDS <- plotDimensions(experiment, assay = "genes", condition = "disease", title = "MDS plot for genes") mirnaMDS <- plotDimensions(experiment, assay = "microRNA", condition = "disease", title = "MDS plot for miRNAs") ggpubr::ggarrange(geneMDS, mirnaMDS, nrow = 1, labels = "AUTO", common.legend = TRUE) ## ----model-------------------------------------------------------------------- ## design the linear model for both genes and miRNAs model <- ~ disease + patient ## ----diffexp------------------------------------------------------------------ ## perform differential expression for genes experiment <- performGeneDE(experiment, group = "disease", contrast = "PTC-NTH", design = model, pCutoff = 0.01) ## perform differential expression for miRNAs experiment <- performMirnaDE(experiment, group = "disease", contrast = "PTC-NTH", design = model, pCutoff = 0.01) ## ----accessDe----------------------------------------------------------------- ## access DE results for genes deGenes <- geneDE(experiment) ## access DE results for miRNAs deMirnas <- mirnaDE(experiment) ## ----volcano, fig.wide=TRUE, fig.cap="Volcano plots of gene and miRNA differential expression. (A) shows the differentially expressed genes, while (B) displays differentially expressed miRNAs."---- ## create a volcano plot for genes geneVolcano <- plotVolcano(experiment, assay = "genes", title = "Gene differential expression") ## create a volcano plot for miRNAs mirnaVolcano <- plotVolcano(experiment, assay = "microRNA", title = "miRNA differential expression") ## plot graphs side by side ggpubr::ggarrange(geneVolcano, mirnaVolcano, nrow = 1, labels = "AUTO", common.legend = TRUE) ## ----thyroidBars, fig.wide=FALSE, fig.cap="Differential expression bar plots for different thyroid genes. Differential expression analysis demonstrated how TG, TPO, DIO2 and PAX8 result downregulated in thyroid cancer."---- ## create a bar plot for all thyroid features plotDE(experiment, features = c("TG", "TPO", "DIO2", "PAX8"), graph = "barplot", linear = FALSE) ## ----species------------------------------------------------------------------ ## list available species for Reactome database supportedOrganisms("Reactome") ## ----ora---------------------------------------------------------------------- ## perform over-representation analysis with GO ora <- enrichGenes(experiment, method = "ORA", database = "GO", organism = "Homo sapiens") ## ----gsea--------------------------------------------------------------------- ## set seed for reproducible results set.seed(1234) ## perform gene-set enrichment analysis with KEGG gse <- enrichGenes(experiment, method = "GSEA", database = "KEGG", organism = "Homo sapiens") ## ----gseaTab, echo=FALSE------------------------------------------------------ ## display the results of GSEA knitr::kable(enrichmentResults(gse), digits = 2, caption = "GSEA results. A table listing the affected KEGG pathways in thyroid cancer.") ## ----oraDot, fig.wide=FALSE, fig.cap="ORA results for downregulated genes. The enrichment of downregulated genes through the gene sets provided by GO database."---- ## create a dot plot for ORA enrichmentDotplot(ora$downregulated, title = "Depleted functions") ## ----gsePlot, fig.wide=FALSE, fig.cap="GSEA-style plot for Thyroid hormone synthesis. This type of plot shows the running sum that GSEA uses to determinate the enrichment score for each pathway."---- ## create a GSEA plot gseaPlot(gse, "Thyroid hormone synthesis", rankingMetric = TRUE) ## ----targets, results='hide', eval=FALSE-------------------------------------- # ## retrieve miRNA target genes # experiment <- getTargets(experiment) ## ----targetsTMP, include=FALSE------------------------------------------------ experiment <- loadExamples() ## ----correlation-------------------------------------------------------------- ## perform a correlation analysis experiment <- mirnaIntegration(experiment, test = "correlation") ## ----correlationResults------------------------------------------------------- ## extract correlation results integrationResults <- integration(experiment) ## take a look at correlation results head(integrationResults) ## ----corPlot, fig.wide=TRUE, fig.cap="Correlation between miRNAs and key thyroid genes. These plots suggest that the upregulation of miR-146b-5p and miR-146b-3p may be responsible for the downregulation of DIO2 and PAX8, respectively."---- ## plot the correlation between miR-146b-5p and DIO2 cor1 <- plotCorrelation(experiment, mirna = "hsa-miR-146b-5p", gene = "DIO2", condition = "disease") ## plot the correlation between miR-146b-3p and PAX8 cor2 <- plotCorrelation(experiment, mirna = "hsa-miR-146b-3p", gene = "PAX8", condition = "disease") ## plot graphs side by side ggpubr::ggarrange(cor1, cor2, nrow = 1, labels = "AUTO", common.legend = TRUE) ## ----association-------------------------------------------------------------- ## perform a one-sided inverse association exp.association <- mirnaIntegration(experiment, test = "association", associationMethod = "fisher-midp", pCutoff = 0.2, pAdjustment = "none") ## ----fry---------------------------------------------------------------------- ## perform the integrative analysis through 'fry' method exp.fry <- mirnaIntegration(experiment, test = "fry", pAdjustment = "none") ## ----sessionInfo, echo=FALSE-------------------------------------------------- sessionInfo()