---
title: "E-value in DNA methylation studies"
output:
  prettydoc::html_pretty:
    theme: cayman
    highlight: github
vignette: >
  %\VignetteIndexEntry{E-value in DNA methylation studies}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```{css, echo=FALSE}
p {
  font-size: 12px;
}

<style type="text/css">
.main-container {
  max-width: 100% !important;
  margin: auto;
}
</style>
```

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# Introductions

In this package, we provide e-value for four DMR (differentially methylated region) detection tools (MethylKit, Metilene, BiSeq and DMRfinder) and general purpose. For `DMR` (`methylKit`, `biseq`, `DMRfinder` or `metilene`), the met-evalue calculation is conducted by the `metevalue.[DMR]` function.

| DMR       | Method                | methyrate Example                   | Method.output Example          |
|:-----------------|:-----------------|:-----------------|:-----------------|
| MethylKit | `metevalue.methylKit` | `data(demo_methylkit_methyrate)`    | `data(demo_methylkit_met_all)` |
| BiSeq     | `metevalue.biseq`     | `data(demo_biseq_methyrate)`        | `data(demo_biseq_DMR)`         |
| DMRfinder | `metevalue.DMRfinder` | `data(demo_DMRfinder_rate_combine)` | `data(demo_DMRfinder_DMRs)`    |
| Metilene  | `metevalue.metilene`  | `data(demo_metilene_input)`         | `data(demo_metilene_out)`      |

Two routines are supported to calculate the combined e-value:

-   Call by **files**: Here the `files` include the outputs of given `DMR` packages and its corresponding e-value of each regions;
-   Call by **R data frames**: Here the `R data frames` are corresponding `data.frame` objects.

Details are as follows.

## Call by files

We design the `metevalue.[DMR]` function to accept similar parameter patterns:

```{r eval=FALSE}
metevalue.[DMR](
  methyrate,                # Output file name of [DMR]
  [DMR].output,             # Output file name of [DMR] with e-value of each region
  adjust.methods = "BH",    # Adjust methods of e-value
  sep = "\t",               # seperator, default is the TAB key.
  bheader = FALSE           # A logical value indicating whether the [DMR].output file
                            # contains the names of the variables as its first line
)
```

Here `[DMR]` could be one of `methylKit`, `biseq`, `DMRfinder` or `metilene`.

## Call by R data frames

We provide the `evalue_buildin_var_fmt_nm` and `varevalue.metilene` function to handle the general DMR e-value calculation in DNA methylation studies:

```{r eval=FALSE}
# Here  `[DMR]` coudle be one of `methylKit`, `biseq`, `DMRfinder` or `metilene`.
method_in_use = "[DMR]"
result = evalue_buildin_var_fmt_nm(
          methyrate,              # Data frame of the methylation rate
          DMR_evalue_output,      # Data frame of output data corresponding to the
                                  # "method" option
          method = method_in_use) # DMR: "metilene", "biseq", "DMRfinder" or "methylKit"
result = list(a = result$a,
              b = result$b,
              a_b = evalue_buildin_sql(result$a, result$b, method = method_in_use))
result = varevalue.metilene(result$a, result$b, result$a_b)
```

Replace `[DMR]` to one of `methylKit`, `biseq`, `DMRfinder` or `metilene` accordingly.

> Notice: for different `[DMR]`, the `data.frame` schemas are **different**!!! Check the R help document for details. Check the [Demo data](#demo-data) section for details. \# Examples

## MethylKit Example

methylKit is a R package for DNA methylation analysis and annotation from high-throughput bisulfite sequencing. The package is designed to deal with sequencing data from RRBS and its variants, but also target-capture methods and whole genome bisulfite sequencing.

Currently, `metevalue` package supports the e-value analytics of the `methylKit` output file.

```{r eval=FALSE}
library(metevalue)

####Simulation Data ####
set.seed(1234)

# methylKit is a Bioconductor package
library(methylKit)
file.list=list( system.file("extdata", 
                            "test1.myCpG.txt", package = "methylKit"),
                system.file("extdata",
                            "test2.myCpG.txt", package = "methylKit"),
                system.file("extdata", 
                            "control1.myCpG.txt", package = "methylKit"),
                system.file("extdata", 
                            "control2.myCpG.txt", package = "methylKit") )


# read the files to a methylRawList object: myobj
myobj=methRead(file.list,
               sample.id=list("test1","test2","ctrl1","ctrl2"),
               assembly="hg18",
               treatment=c(1,1,0,0),
               context="CpG"
)

meth=unite(myobj, destrand=FALSE)
meth.C <- getData(meth)[,seq(6,ncol(meth),3)]
meth.T <- getData(meth)[,seq(7,ncol(meth),3)]
mr <- meth.C/(meth.C + meth.T)
chr_pos = getData(meth)[,1:2]
methyrate = data.frame(chr_pos,mr)
names(methyrate) = c('chr', 'pos', rep('g1',2), rep('g2',2))
region<-tileMethylCounts(myobj)
meth<-unite(region,destrand=F)
myDiff<-calculateDiffMeth(meth)
met_all<-getMethylDiff(myDiff,type="all")

example_tempfiles = tempfile(c("rate_combine", "methylKit_DMR_raw"))
tempdir()
write.table(methyrate, file=example_tempfiles[1], row.names=F, col.names=T, quote=F, sep='\t')
write.table (met_all, file=example_tempfiles[2], sep ="\t", row.names =F, col.names =T, quote =F)
```

`evalue.methylKit` function could be used to tackle the problem.

```{r eval=FALSE}
result = metevalue.methylKit(example_tempfiles[1], example_tempfiles[2], bheader = T)
str(result)
```

Alternatively, one could use the build-in functions to derive functions which avoid the file operation:

```{r eval=FALSE}
result = evalue_buildin_var_fmt_nm(methyrate, met_all, method="methylKit")
result = list(a = result$a, 
              b = result$b, 
              a_b = evalue_buildin_sql(result$a, result$b, method="methylKit"))
result = varevalue.metilene(result$a, result$b, result$a_b)
str(result)
```

## BiSeq Example

First, we load the methylation data at CpG site levels from 'BiSeq' package. Then we clustered CpG sites into DMRs using 'BiSeq'.

```{r eval=FALSE}
library(BiSeq)
library(dplyr)
data(rrbs)
rrbs.rel <- rawToRel(rrbs)
methyrate <- methLevel(rrbs.rel)
methyrate <- data.frame(methyrate)
methyrateq = cbind(rows = as.numeric(row.names(methyrate)), methyrate)
methypos = data.frame(rows = as.numeric(row.names(methyrate)), rowRanges(rrbs))
methyrate = left_join(methypos, methyrateq)
methyrate = methyrate[,c(2,3,7:16)]
names(methyrate) <- c('chr','pos',rep('g1',5),rep('g2',5))

rrbs.clust.unlim <- clusterSites(object = rrbs,perc.samples = 3/4,min.sites = 20,max.dist = 100)

clusterSitesToGR(rrbs.clust.unlim)
ind.cov <- totalReads(rrbs.clust.unlim) > 0

quant <- quantile(totalReads(rrbs.clust.unlim)[ind.cov])
rrbs.clust.lim <- limitCov(rrbs.clust.unlim, maxCov = quant)
predictedMeth <- predictMeth(object = rrbs.clust.lim)

test<- predictedMeth[, colData(predictedMeth)$group == "test"]
control <- predictedMeth[, colData(predictedMeth)$group == "control"]
mean.test <- rowMeans(methLevel(test))
mean.control <- rowMeans(methLevel(control))

betaResults <- betaRegression(formula = ~group,link = "probit",object = predictedMeth,type = "BR")
vario <- makeVariogram(betaResults)
vario.sm <- smoothVariogram(vario, sill = 0.9)

locCor <- estLocCor(vario.sm)
clusters.rej <- testClusters(locCor)
clusters.trimmed <- trimClusters(clusters.rej)
DMRs <- findDMRs(clusters.trimmed,max.dist = 100,diff.dir = TRUE)


example_tempfiles = tempfile(c('rate_combine', 'BiSeq_DMR'))
write.table(methyrate, example_tempfiles[1], row.names=F, col.names=T, quote=F, sep='\t')
write.table(DMRs, example_tempfiles[2], quote=F, row.names = F,col.names = F, sep = '\t')
```

Finally, we added E-values and adjusted E-values as additional columns to the output file of 'BiSeq'.`BiSeq_evalue` function could be used to tackle the problem.

```{r eval=FALSE}
hh <- data.frame(DMRs)
result = evalue_buildin_var_fmt_nm(rate_combine, hh, method="biseq")
result = list(a = result$a,  b = result$b, 
              a_b = evalue_buildin_sql(result$a, result$b,method="biseq"))
result = varevalue.metilene(result$a, result$b, result$a_b)
str(result)
```

## DMRfinder Example

Given the input file

-   `rate_combine_DMRfinder`: a file containing methylation rates at each CpG site

-   `DMRfinder_DMR`: the output file from 'DMRfinder'

```{r eval=FALSE}
rate_combine <- read.table("rate_combine_DMRfinder", header = T)
head(rate_combine)

DMRs <- read.table("DMRfinder_DMR", header = T)
head(DMRs)

```

Adding E-values and adjusted E-values as additional columns to file 'DMRfinder_DMR'

```{r eval=FALSE}
result <- evalue.DMRfinder('rate_combine_DMRfinder', 'DMRfinder_DMR')
head(result)
```

Alternatively, function `varevalue.metilene` can also provides e-value and adjusted e-value.

```{r eval=FALSE}
result = evalue_buildin_var_fmt_nm(rate_combine, DMRs, method="DMRfinder")
result = list(a = result$a, 
              b = result$b, 
              a_b = evalue_buildin_sql(result$a, result$b, method="DMRfinder"))
result = varevalue.metilene(result$a, result$b, result$a_b)
head(result)
```

## Metilene Example

Given

-   `metilene.input`: the input file of `Metilene` containing methylation rates at each CpG site
-   `metilene.out`: the output file of `Metilene`

```{r eval=FALSE}
input <- read.table("metilene.input", header = T)
head(input)

out <- read.table("metilene.out", header = F)
head(out)

```

Adding E-values and adjusted E-values as additional columns to `metilene.out`

```{r eval=FALSE}
result <- evalue.metilene('metilene.input', 'metilene.out')
head(result)
```

Alternatively, function `varevalue.metilene` can also provides e-value and adjusted e-value.

```{r eval=FALSE}
result = evalue_buildin_var_fmt_nm(input, out, method="metilene")
result = list(a = result$a, 
              b = result$b, 
              a_b = evalue_buildin_sql(result$a, result$b, method="metilene"))
result = varevalue.metilene(result$a, result$b, result$a_b)
head(result)
```

# Misc

## Demo data {#demo-data}

Demo data for different `metevalue.[DMR]` functions are listed in the section.

### Input Data Examples: MethylKit

**methyrate Example**

|chr   |     pos|        g1|        g1|        g2|        g2|
|:-----|-------:|---------:|---------:|---------:|---------:|
|chr21 | 9853296| 0.5882353| 0.8048048| 0.8888889| 0.8632911|
|chr21 | 9853326| 0.7058824| 0.7591463| 0.8750000| 0.7493404|


- chr: string Chromosome
- pos: int Position
- g1~g2: methylation rate data in groups, repeat 8 times. Notice that there are "NaN" within the feature columns.


**methylKit.output Example**

|chr| start|     end| strand | pvalue| qvalue|meth.diff|
|:-----|-------:|-------:|:------|------:|------:|---------:|
|chr21 | 9927001| 9928000|*      |      0|      0| -34.07557|
|chr21 | 9944001| 9945000|*      |      0|      0| -40.19089|


- chr: Chromosome
- start: The positions of the start sites of the corresponding region
- end: The positions of the end sites of the corresponding region
- strand: Strand
- pvalue: The adjusted p-value based on BH method in MWU-test
- qvalue: cutoff for qvalue of differential methylation statistic
- methyl.diff: The difference between the group means of methylation level


### Input Data Examples: BiSeq

**methyrate Example**

|chr  |    pos|        g1| g1|  g1|  g1|  g1|     g2|        g2|        g2|        g2|        g2|
|:----|------:|---------:|--:|---:|---:|---:|------:|---------:|---------:|---------:|---------:|
|chr1 | 870425| 0.8205128|  1| 0.7| NaN| NaN| 0.3125| 0.7419355| 0.2461538| 0.1794872| 0.2413793|
|chr1 | 870443| 0.8461538|  1| 0.7| NaN| NaN| 0.3750| 0.3225806| 0.2923077| 0.0512821| 0.2413793|


- chr: Chromosome
- pos: int Position
- g1~g2: methylation rate data in groups


**biseq.output Example**

|seqnames |  start|    end| width|strand |  median.p| median.meth.group1| median.meth.group2| median.meth.diff|
|:--------|------:|------:|-----:|:------|---------:|------------------:|------------------:|----------------:|
|chr1     | 872369| 872616|   248|*      | 0.0753559|          0.9385462|          0.8666990|        0.0710524|
|chr1     | 875227| 875470|   244|*      | 0.0000026|          0.5136315|          0.1991452|        0.2942668|


- seqnames: Chromosome
- start: The positions of the start sites of the corresponding region
- end: The positions of the end sites of the corresponding region
- width: The range of the corresponding region
- strand: Strand
- median.p: The median of p-values in the corresponding region
- median.meth.group1 : The median of methylation level for the corresponding segment of group 1
- median.meth.group2 : The median of methylation level for the corresponding segment of group 2
- median.meth.diff: The median of the difference between the methylation level



### Input Data Examples: DMRfinder

**methyrate Example**

|chr  |       pos| g1|      g1.1| g2| g2.1|
|:----|---------:|--:|---------:|--:|----:|
|chr1 | 202833315|  0| 0.0000000|  0|    0|
|chr1 | 202833323|  1| 0.8095238|  1|    1|


- chr: Chromosome
- pos: int Position
- g1~g2: methylation rate data in groups

**DMRfinder.output Example**

|chr   |    start|      end| CpG| Control.mu|  Exptl.mu| Control.Exptl.diff| Control.Exptl.pval|
|:-----|--------:|--------:|---:|----------:|---------:|-------------------:|-------------------:|
|chr8  | 25164078| 25164102|   3|  0.9241646| 0.7803819|          -0.1437827|           0.0333849|
|chr21 |  9437432|  9437538|  14|  0.7216685| 0.1215506|          -0.6001179|           0.0000000|


- chr: Chromosome
- start: The positions of the start sites of the corresponding region
- end: The positions of the end sites of the corresponding region
- CpG: The number of CpG sites within the corresponding region
- Control.mu: The average methylation rate in control group
- Expt1.mu: The average methylation rate in experiment group
- Control.Expt1.diff: The methylation difference between control and experiment groups
- Control.Expt1.pval: P-value based on Wald-test


### Input Data Examples: DMRfinder

**methyrate Example**

|chrom |     pos|        g1| g1.1|      g1.2| g1.3| g1.4|      g1.5| g1.6|      g1.7|        g2| g2.1| g2.2| g2.3|      g2.4| g2.5| g2.6| g2.7|
|:-----|-------:|---------:|----:|---------:|----:|----:|---------:|----:|---------:|---------:|----:|----:|----:|---------:|----:|----:|----:|
|chr21 | 9437433| 0.9285714|   NA| 0.7222222| 0.75|    1| 0.6666667|    1| 0.8695652| 0.0000000|    0|    0|    0| 0.0000000|  0.0|   NA| 0.00|
|chr21 | 9437445| 1.0000000|   NA| 0.9444444| 0.75|    1| 0.6666667|    0| 0.8695652| 0.6111111|    0|    0|    0| 0.7333333|  0.6|   NA| 0.75|


- chr: Chromosome
- pos: int Position
- g1~g2: methylation rate data in groups


**metilene.output Example**

| chr  | start  |  end  |  q-value | methyl.diff  |  CpGs  |  p | p2  | m1 | m2 |
|:-----|-------:|-------:|--:|--------:|--:|--:|--:|-------:|-------:|
|chr21 | 9437432| 9437540|  0| 0.610989| 26|  0|  0| 0.73705| 0.12606|
|chr21 | 9708982| 9709189|  0| 0.475630| 28|  0|  0| 0.58862| 0.11299|

- chr: Chromosome
- start: The positions of the start sites of the corresponding region
- end: The positions of the end sites of the corresponding region
- q-value: The adjusted p-value based on BH method in MWU-test
- methyl.diff: The difference between the group means of methylation level
- CpGs: The number of CpG sites within the corresponding region
- p : p-value based on MWU-test
- p2: p-value based on 2D KS-test
- m1: The absolute mean methylation level for the corresponding segment of group 1
- m2: The absolute mean methylation level for the corresponding segment of group 2
