Adding updated documentation for mlml

Author: andrewdavidsmith

docs/content/mlml.md

@@ -1,109 +1,168 @@
 # mlml - estimation of hydroxymethylation and methylation levels
 
 ## Synopsis
-```shell
-$ dnmtools mlml [OPTIONS] -u <input-bs.meth> -m <input-oxbs.meth> -h <input-tab.meth>
+```console
+$ dnmtools mlml [OPTIONS] -u <input-bs-seq> -m <input-oxbs-seq> -h <input-tab-seq>
 ```
 
 ## Description
 
 If you are interested in estimating hydroxymethylation level and have
 any two of Tet-Assisted Bisulfite sequencing (TAB-seq), oxidative
-bisulfite sequencing (oxBS-seq) and BS-seq data available, you can use
-`mlml`  to perform consistent and simultaneous estimation.
-
-The input file format could be the default
-[counts](../counts) output format,
-or BED format file with 6 columns as the example below:
+bisulfite sequencing (oxBS-seq) or ordinary bisulfite sequencing
+(BS-seq) data available, `mlml` can perform simultaneous estimation of
+5mC and 5hmC levels, separately for each site.
 
+The input file format is the "counts" format, which is generated by
+the dnmtools [counts](../counts) command after mapping reads. This
+is the general appearance of each line in the input files:
 ```txt
-chr1  3001345  3001346  CpG:9  0.777777777778  +
+chr1    3001345    +    CpG    0.777777777778    9
 ```
-
-Here the fourth column indicates that this site is a CpG site, and the
-number of reads covering this site is 9. The fifth column is the
-methylation level of the CpG site, ranging from 0 to 1. Note that all
-input files must be sorted. Assume you have three input files ready:
-`input-bs-seq.meth`, `input-oxbs-seq.meth` and `input-tab-seq.meth`.
-The following command will take all the inputs:
-
-```shell
-$ dnmtools mlml -v -u input-bs-seq.meth -m input-oxbs-seq.meth -h input-tab-seq.meth -o result.txt
+In general `mlml` works on CpG sites. There is no reason you cannot
+use it on other sites, but we've only ever tested it on CpG sites.
+
+All input files for `mlml` must be sorted the same way, and have
+corresponding sites -- each must have the same number of sites, at the
+same positions within the genome, and in the same order.
+
+Assume you have three input files: `input-bs-seq.counts`,
+`input-oxbs-seq.counts` and `input-tab-seq.counts`.  Each of these
+files would be produced by ordinary BS-seq, oxBS-seq and TAB-seq,
+respectively.  The following command will use these inputs:
+```console
+$ dnmtools mlml -u input-bs-seq.counts -m input-oxbs-seq.counts -h input-tab-seq.counts
 ```
+The results for the above command will be printed to the terminal with
+one line per site (line) in the input files -- that might be a
+lot. One line will look something like this:
+```txt
+chr1    300    +     CpG    0.635788    0.259378    0.104835    0
+```
+Each line has 8 columns, the first 4 of which correspond to those in
+the [counts](../counts) format indicated above. These are
+*tab-delimited*.  Here is the meaning of the rest of the columns:
+
+1. chromosome
+2. position
+3. strand
+4. label
+5. 5mC: estimated fraction of molecules with the 5mC mark
+6. 5hmC: estimated fraction of molecules with the 5mC mark
+7. neither: estimated fraction with neither 5mC nor 5hmC
+8. number of conflicting sites: this integer value is explained below
+
+If only two types of input are available, e.g. `input-bs-seq.counts`
+and `input-oxbs-seq.counts`, then the following command would be
+appropriate:
+```console
+$ dnmtools mlml -u input-bs-seq.counts -m input-oxbs-seq-seq.counts -o result.txt
+```
+The output would be written to the file `results.txt` because of the
+`-o` argument. Obviously more data is better, and estimates from `mlml`
+are more accurate with the 3 types of data. But in the vast majority
+of cases only one of oxBS-seq or TAB-seq will be available. The
+arguments `-u`, `-m` and `-b` are used to indicate the following:
 
-If only two types of input are available, e.g. `input-bs-seq.meth and
-`input-oxbs-seq.meth`, then use the following command:
+* -u: BS-seq input, think of it as directly informing about "unmethylated"
+* -m: oxBS-seq input, "m" because it gives exactly 5mC information.
+* -h: TAB-seq input, "h" because it gives exactly 5hmC information.
 
-```shell
-$ dnmtools mlml -u input-bs-seq.meth -m input-oxbs-seq-seq.meth -o result.txt
-```
+You must provide at least 2 of these, and be sure to use the
+appropriate argument for each data type. `mlml` works by taking
+advantage of the properties of each data type, so swapping those will
+result in useless output.
 
 In some cases, you might want to specify the convergence tolerance for
 EM algorithm. This can be done through `-t` option.  For example:
-
-```shell
-$ dnmtools mlml -u input-bs-seq.meth -m input-oxbs-seq.meth -o result.txt -t 1e-2
+```console
+$ dnmtools mlml -u bs-seq.counts -m oxbs-seq.counts -o result.txt -t 1e-2
 ```
-
 This command will make the iteration process stop when the difference
-of estimation between two iterations is less than 10^−2 . Thea value
-format can be scientific notation, e.g. 1e-5, or float number, e.g.
-0.00001.
-
-The output of `mlml` is tab-delimited format. Here is an example:
-
-```txt
-chr11   15      16      0.166667        0.19697 0.636364        0
-chr12   11      12      0.222222        0       0.777778        2
+of estimation between two iterations is less than 1e-02. The value
+format can be scientific notation, e.g. 1e-5, or as a fixed point
+fraction, e.g.  0.00001.
+
+Conflicts: To calculate the last column in the output, a binomial test
+is performed for each input methylation level (either 2 or 3 in total
+depending on the number of input files). The tests are done separately
+for estimated maximum likelihood level (for each of 5mC, 5hmC and
+"neither"). In each case, check if the estimate falls outside the
+confidence interval calculated using coverage and observed fraction
+from the input file. If more than one estimate at a given site falls
+outside the confidence interval, we flag it as conflicting. This will
+happen by random chance, but if the same sites show conflicts in
+multiple independent analyses, it might indicate that the sites
+themselves are unreliable in one of the types of experiments.
+
+## Examples
+
+Assume the data files are as follows:
+```console
+$ cat bs.counts
+chr1    100     +       CpG     0.95    20
+chr1    200     +       CpG     0.8     20
+chr1    300     +       CpG     0.9     30
+chr1    400     +       CpG     0.5     2
+$ cat oxbs.counts
+chr1    100     +       CpG     0.8    20
+chr1    200     +       CpG     0.6    20
+chr1    300     +       CpG     0.6    10
+chr1    400     +       CpG     0.6    10
+$ cat tab.counts
+chr1    100     +       CpG     0.15    20
+chr1    200     +       CpG     0.2     30
+chr1    300     +       CpG     0.2     5
+chr1    400     +       CpG     0.1     20
 ```
-
-The columns are chromosome name, start position, end position, 5-mC
-level, 5-hmC level, unmethylated level and number of conflicts. To
-calculate the last column, a binomial test is performed for each input
-methylation level (can be 2 or 3 in total depending on parameters). If
-the estimated methylation level falls out of the confidence interval
-calculated from input coverage and methylation level, then such event
-is counted as one conflict. It is recommended to filter estimation
-results based on the number of conflicts; if more conflicts happens on
-one site then it is possible that information from such site is not
-reliable.
+Then running `mlml` would give the following results:
+```console
+$ dnmtools mlml -bsseq bs.counts -tabseq tab.counts -oxbsseq oxbs.counts
+chr1    100 +   CpG 0.8         0.15        0.05        0
+chr1    200 +   CpG 0.6         0.2         0.2         0
+chr1    300 +   CpG 0.635788    0.259378    0.104835    0
+chr1    400 +   CpG 0.565183    0.0939689   0.340848    0
+```
+You should be able to calculate the estimates for the 1st and 2nd
+output lines easily in your head. The third and fourth lines are more
+difficult.
 
 ## Options
 
 ```txt
  -o, -output
 ```
-output file (default: STDOUT)
+output file (default: standard output)
 
 ```txt
  -u, -bsseq
 ```
-input BS-seq methcounts file
+input file for BS-seq counts
 ```txt
  -h, -tabseq
 ```
-input TAB-seq methcounts file
+input file for TAB-seq counts
 ```txt
  -m, -oxbsseq
 ```
-input oxBS-seq methcounts file
+input file for oxBS-seq counts
 ```txt
  -t, -tolerance
 ```
-EM convergence threshold (default: 0.000000)
+convergence tolerance (default: 1e-10)
 ```txt
  -a, -alpha
 ```
-significance level of binomial test for each site (default: 0.050000)
+significance level of binomial test at each site (default: 0.05)
 ```txt
  -H, -outh
 ```
-hmC pseudo methcount output file (default: null)
+5hmC estimated levels at each site (default: none)
 ```txt
  -M, -outm
 ```
-mC pseudo methcount output file (default: null)
+5mC estimated levels at each site (default: none)
 ```txt
  -v, -verbose
 ```
-print more run info to STDERR while the program is running.
+print more information while running