This document guides one through all available functions of the anamiR
package. Package anamiR aims to find potential miRNA-target gene interactions from both miRNA and mRNA expression data.
Traditional miRNA analysis method is to use online databases to predict miRNA-target gene interactions. However, the inconsistent results make interactions are less reliable. To address this issue, anamiR integrates the whole expression analysis with expression data into workflow, including normalization, differential expression, correlation and then databases intersection, to find more reliable interactions. Moreover, users can identify interactions from interested pathways or gene sets.
anamiR is on Bioconductor and can be installed following standard installation procedure.
install.packages("BiocManager")
BiocManager::install("anamiR")
To use,
library(anamiR)
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
## Registered S3 method overwritten by 'openssl':
## method from
## print.bytes Rcpp
## Setting options('download.file.method.GEOquery'='auto')
## Setting options('GEOquery.inmemory.gpl'=FALSE)
The general workflow can be summarized as follows,
Basically there are six steps, corresponding to six R functions, to complete the whole analysis:
As shown in the workflow, not only samples of paired miRNA and mRNA expression data, but phenotypical information of miRNA and mRNA are also required for the analysis. Since anamiR reads data in expression matrices, data sources are platform and technology agnostic. particularly, expression data from microarray or next generation sequencing are all acceptable for anamiR. However, this also indicates the raw data should be pre-processd to the expression matrices before using anamiR.
Columns for samples. Rows for genes
GENE SmapleA SamplB ...
A 0.1 0.2
B -0.5 -0.3
C 0.4 0.1
Columns for samples. Rows for miRNAs
miRNA SampleA SampleB ...
A 0.1 0.5
B -0.5 2.1
C 0.4 0.3
Cloumns for samples. Rows for feature name, including two groups, multiple groups, continuous data.
Feature groupA groupA groupA ...
SampleA 123.33 A A
SampleB 120.34 B C
SampleC 121.22 A B
Now, we show an example using internal data for anamiR workflow.
To demonstrate the usage of the anamiR
package, the package contains 30 paired
miRNA and mRNA breast cancer samples, which are selected from 101 miRNA samples and
114 mRNA samples from GSE19536. As for phenotype data (hybridization information),
there are three types of information in it, including two-groups, multi-groups,
continuous data.
The mRNA data was conducted by Agilent-014850 Whole Human Genome Microarray 4x44K G4112F (Probe Name version) platform and the miRNA data was generated from Agilent-019118 Human miRNA Microarray 2.0 G4470B (miRNA ID version).
First of all, load the internal data and check the format.
data(mrna)
data(mirna)
data(pheno.mirna)
data(pheno.mrna)
Basically, the format of input data should be the same as the internal data. For mRNA expression data should be like,
mrna[1:5, 1:5]
## BC.M.014 BC.M.015 BC.M.017 BC.M.019 BC.M.023
## TRIM35 8.725157 8.699302 9.864870 8.986249 9.253038
## NR2C2AP 10.973569 11.318088 11.212974 11.803147 11.823159
## CENPO 7.404785 8.417221 8.272822 8.670029 9.392538
## ADCY8 5.785176 6.066669 6.240934 6.345486 6.037214
## RGS22 8.828413 7.844503 7.712051 7.078598 8.712837
As for miRNA expression data,
mirna[1:5, 1:5]
## BC.M.014 BC.M.015 BC.M.017 BC.M.019 BC.M.023
## hsa-miR-139-5p 3.301236 4.072743 3.476371 3.957260 6.662144
## hsa-miR-575 9.525705 10.097954 10.790396 9.365697 9.289803
## hsa-miR-33b* 4.089836 4.297618 4.051852 3.842798 3.044545
## hsa-miR-335 4.146060 5.615910 6.274007 5.655296 9.309553
## hsa-miR-16-2* 3.407055 4.023072 3.052135 4.208022 4.572399
And the phenotype data, (NOTICE:users should arrange the case columns front of the control columns.)
pheno.mrna[1:3, 1:3]
## Subtype Survival ER
## BC.M.014 "breast cancer subtype: Lum A" "126.611842105263" "case"
## BC.M.015 "breast cancer subtype: Lum A" "126.513157894737" "case"
## BC.M.017 "breast cancer subtype: Lum A" "80.9868421052632" "case"
pheno.mrna[28:30, 1:3]
## Subtype Survival ER
## BC.M.221 "breast cancer subtype: Basal" "115.197368421053" "control"
## BC.M.318 "breast cancer subtype: ERBB2" "110.526315789474" "control"
## BC.M.709 "breast cancer subtype: Basal" "15.2960526315789" "control"
Actually, the phenotype data of miRNA and mRNA share the same contents,but in this case, we still make it in two data to prevent users from being confused about it.
Secondly, we normalize data. (If you use the normalized data, you can skip this step.)
se <- normalization(data = mirna, method = "quantile")
For this function, there are three methods provided, including quantile
,
rank.invariant
, normal
. For more detail, Please refer to their
documentation.
Note that internal data have already been normalized, here is only for demonstration.
Before entering the main workflow, we should put our data and phenotype
information into SummarizedExperiment
class first, which you can get
more information from SummarizedExperiment.
mrna_se <- SummarizedExperiment::SummarizedExperiment(
assays = S4Vectors::SimpleList(counts=mrna),
colData = pheno.mrna)
mirna_se <- SummarizedExperiment::SummarizedExperiment(
assays = S4Vectors::SimpleList(counts=mirna),
colData = pheno.mirna)
Third, we will find differential expression genes and miRNAs.
There are three statitical methods in this function. here, we use
t.test
for demonstration.
mrna_d <- differExp_discrete(se = mrna_se,
class = "ER", method = "t.test",
t_test.var = FALSE, log2 = FALSE,
p_value.cutoff = 0.05, logratio = 0.5
)
mirna_d <- differExp_discrete(se = mirna_se,
class = "ER", method = "t.test",
t_test.var = FALSE, log2 = FALSE,
p_value.cutoff = 0.05, logratio = 0.5
)
This function will delete genes and miRNAs (rows), which do not differential express, and add another three columns represent fold change (log2), p-value, adjusted p-value.
Take differential expressed mRNA data for instance,
nc <- ncol(mrna_d)
mrna_d[1:5, (nc-4):nc]
## log-ratio P-Value P-adjust mean_case mean_control
## NR2C2AP -1.0712486 2.316469e-04 5.312839e-04 11.556393 12.627641
## CENPO -0.8558411 2.682517e-04 5.748560e-04 8.349869 9.205710
## RGS22 1.6894848 2.039121e-05 1.631297e-04 8.398516 6.709031
## KIF3C -0.5848147 4.541900e-05 2.270950e-04 8.279248 8.864063
## SENP5 -0.7626540 8.630361e-07 1.917858e-05 10.737840 11.500494
Before using collected databases for intersection with potential miRNA-target gene interactions, we have to make sure all miRNA are in the latest annotation version (miRBase 21). If not, we could use this function to do it.
mirna_21 <- miR_converter(data = mirna_d, remove_old = TRUE,
original_version = 17, latest_version = 21)
Now, we can compare these two data,
# Before
head(row.names(mirna_d))
## [1] "hsa-miR-575" "hsa-miR-150*" "hsa-miR-198" "hsa-miR-342-3p"
## [5] "hsa-miR-342-5p" "hsa-miR-134"
# After
head(row.names(mirna_21))
## [1] "hsa-miR-575" "hsa-miR-150-3p" "hsa-miR-198" "hsa-miR-342-3p"
## [5] "hsa-miR-342-5p" "hsa-miR-134-5p"
Note that user must put the right original version into parameter, because it is an important information for function to convert annotation.
To find potential miRNA-target gene interactions, we should
combine the information in two differential expressed data,
which we obtained from differExp_discrete
.
cor <- negative_cor(mrna_data = mrna_d, mirna_data = mirna_21,
method = "pearson", cut.off = -0.5)
For output,
head(cor)
## miRNA Gene Correlation( -0.5 ) logratio(miRNA)
## [1,] "hsa-miR-198" "NR2C2AP" "-0.60749521052608" "1.434328524"
## [2,] "hsa-miR-198" "SENP5" "-0.54646311346452" "1.434328524"
## [3,] "hsa-miR-134-5p" "SENP5" "-0.508211941696094" "1.0180378"
## [4,] "hsa-miR-198" "AKAP1" "-0.540628213927599" "1.434328524"
## [5,] "hsa-miR-342-3p" "PIF1" "-0.565476513859015" "2.19495706"
## [6,] "hsa-miR-342-5p" "PIF1" "-0.621317621065092" "2.115864648"
## P-adjust(miRNA) mean_case(miRNA) mean_control(miRNA)
## [1,] "0.0287086140383543" "4.963163064" "3.52883454"
## [2,] "0.0287086140383543" "4.963163064" "3.52883454"
## [3,] "0.000137286128508541" "8.14071616" "7.12267836"
## [4,] "0.0287086140383543" "4.963163064" "3.52883454"
## [5,] "1.07192166300769e-09" "11.34977766" "9.1548206"
## [6,] "9.42861994530635e-11" "7.135666308" "5.01980166"
## logratio(gene) P-adjust(gene) mean_case(gene)
## [1,] "-1.071248568" "0.000531283948592135" "11.556392652"
## [2,] "-0.762654017199999" "1.91785808708134e-05" "10.7378403228"
## [3,] "-0.762654017199999" "1.91785808708134e-05" "10.7378403228"
## [4,] "-1.2822530884" "0.000659175618505387" "12.1560164316"
## [5,] "-1.5501717548" "0.00020239427080686" "8.5069645972"
## [6,] "-1.5501717548" "0.00020239427080686" "8.5069645972"
## mean_control(gene)
## [1,] "12.62764122"
## [2,] "11.50049434"
## [3,] "11.50049434"
## [4,] "13.43826952"
## [5,] "10.057136352"
## [6,] "10.057136352"
As the showing list
, each row is a potential interaction,
and only the row that correlation coefficient < cut.off would
be kept in list.
Note that in our assumption, miRNAs negatively regulate expression of their target genes, that is, cut.off basically should be negative decimal.
There is a function for user to see the heatmaps about the miRNA-target gene interactions remaining in the correlation analysis table.
heat_vis(cor, mrna_d, mirna_21)