## ----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)

## ----seach_disease------------------------------------------------------------
## identify the EFO ID corresponding to antidepressant response
searchDisease("antidepressant")

## ----snp_association----------------------------------------------------------
## detect SNPs occuring at DE-miRNAs loci
association <- findMirnaSNPs(experiment, "response to antidepressant")

## ----trackPlot, fig.dim=c(7, 3.5), fig.cap="Track plot for miRNA-SNPs. This trackplot shows the proximity of rs2402960 with the locus that encodes for miR-182."----
## create a track plot to represent disease-SNPs at DE-miRNA loci
mirVariantPlot("rs2402960", association, showSequence = FALSE)

## ----snpEvidence--------------------------------------------------------------
## retrieve the evidence supporting SNP-trait association
snpEvidence <- getEvidence(variant = "rs2402960",
                           diseaseEFO = "response to antidepressant")

## take a look at the evidence table
head(snpEvidence)

## ----targets, results='hide'--------------------------------------------------
## retrieve miRNA target genes
experiment <- getTargets(experiment)

## ----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")

## ----intTargEnr, fig.cap="Functional enrichment of integrated targets. This dot plot shows the enriched diseases for downregulated genes."----
## enrichment of integrated targets
oraTarg <- enrichTargets(experiment, database = "DO")

## show a dot plot for the enrichment of downregulated genes
enrichmentDotplot(oraTarg$downregulated,
                  showTerms = 7,
                  showTermsParam = "padj",
                  title = "Depleted diseases")

## ----augmented_networks-------------------------------------------------------
## create miRNA-augmented networks using KEGG pathways
networks <- preparePathways(experiment,
                            database = "KEGG",
                            organism = "Homo sapiens",
                            minPc = 20)

## ----taipa--------------------------------------------------------------------
## only consider a smaller set of augmented networks
networks <- networks[seq(15, 30)]

## set seed for reproducible results
set.seed(1234)

## perform TAIPA with 1000 permutations
taipa <- topologicalAnalysis(experiment,
                             pathways = networks,
                             nPerm = 1000)

## ----integratedPathways-------------------------------------------------------
## extract the results of TAIPA
perturbedNetworks <- integratedPathways(taipa)

## ----integrationDot, fig.wide=FALSE, fig.cap="The perturbation of miRNA-mRNA networks in thyroid cancer. This dot plot display the impairment of thyroid hormone production."----
## produce a dotplot that shows the most affected networks
integrationDotplot(taipa)

## ----thyroidNetwork, fig.wide=TRUE, fig.cap="Impaired network involved in thyroid hormone synthesis. The network created by MIRit suggests that the upregulation of miR-146b-5p and miR-146b-3p may be responsible for diminished expression of PAX8, which in turn causes reduced transcription of thyroid hormone."----
## plot the impaired network responsible for reduced TG synthesis
visualizeNetwork(taipa, "Thyroid hormone synthesis")

## ----sessionInfo, echo=FALSE--------------------------------------------------
sessionInfo()